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1963  A.  Bjerhammar  solved  the  geodetic  boundary  value  problem  by 
applying  Poisson's  integral  equation  for  a finite  set  of  observed  free-air  gravity 
anomalies.  Due  to  the  relation  between  the  number  of  observations  (m)  and  the 
number  of  chosen  unknowns  (N)  different  solutions  are  obtained:  non-singular 
(m  = N),  least  squares  (m  > N)  and  minimum  norm  solutions  (m  < N).  In  the 
special  case  N(^  it  is  shown  that  the  Bjerhammar  solution  with  Poisson's 
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/'’kernel  and  a solution  by  collocation  with  the  corresponding  kernel  are  identical. 
Bjerhammar's  method  is  generalized  by  using  other  kernel  functions,  and  each 
minimum  norm  solution  is  shown  to  correspond  to  one  specific  set  of  degree 
variances  in  collocation. 

The  impulse  approaches  (reflexive  prediction,  Dirac  method)  of 
Bjerhammar  are  presented.  \ The  kernel  function  of  the  non-singular  Dirac 
, method  is  obtained  from  that  in  collocation  by  substituting  c„  /(2n+l)  (rBa/r3r1)1'  + 3 
by  J cn/  (2n+i)  (rB  / r3y  + a.  Hence  the  Dirac  kernel  is  not  symmetric  in  contrast 
to  the  one  in  collocation.  Moreover,  it  is  shown  that  for  a given  radius  ra  of 
the  Bjerhammar  sphere,  the  Dirac  method  gives  a better  conditioning  of  the 
equation  system.  It  is  demonstrated  that  the  two  solutions  are  identic?’,  for 
Poisson's  kernel,  if  the  depth  of  the  Bjerhammar  sphere  in  coll  cation  is  half 
of  that  in  the  application  of  the  Dirac  approach.  The  solutio.  by  collocation  is 
therefore  twice  as  sensitive  to  the  choice  of  radius  as  is  the  latter  method. 

-"In  the  theoretical  case  with  a continuous  coverage  of  observations  at 
the  surface  of  the  earth,  it  is  shown  that  both  the  Dirac  method  and  collocation 
give  a unique  solution  for  any  choice  of  positive  degree  variances  of  the  kernel 
functions,  whenever  the  solutions  exist.  However,  the  intermediate  solutions 
for  Ag*  and  X at  the  Bjerhammar  sphere  do  not  exist  in  general,  if  collocation 
is  applied  by  solving  the  Wiener-Hopf  integral  equation,  a convergent  solution 
is  proved  outside  a sphere.  However,  inside  the  bounding  sphere  of  the  earth 
the  convergence  is  still  not  proved. 
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1.  Introduction 

In  1963  A.  Bjerhammar  formulated  the  boundary  value  problem  of 
physical  geodesy  in  the  following  way:  "A  finite  number  of  gravity  data  (gravity 
anomalies)  is  given  for  a non-spherical  surface  and  it  is  required  to  find  such 
a solution+that  the  boundary  values  are  satisfied  in  all  given  points".  (Discrete 
boundary  value  problem.) 

For  the  solution  of  this  problem  Bjerhammar  used  a spherical  reference 
surface  completely  embedded  in  the  earth  with  its  center  at  the  earth's  center  of 
gravity  and  rotating  with  the  same  angular  velocity  as  the  real  earth.  The 
observed  gravity  anomalies  ( Ag  j)  are  reduced  to  the  internal  sphere  (the 
Bjerhammar  sphere)  by  means  of  Poisson's  integral  equation  for  the  harmonic 
function  rAg  (see  Figure  1): 

P. 


Figure  1.  The  observation  Ag  at  the  surface  of  the  earth  is  reduced  to  the  fictitious 
field  Ag*  on  the  internal  sphere  by  means  of  Poisson's  integral  equation. 


t for  the  height  anomalies  (see  section  5) 


where 


Agj  = real  gravity  - theoretical  gravity 

Ag*  = (fictitious)  gravity  anomaly  at  the  Bjerhammar  sphere 
rj  = geocentric  distance  to  the  point  Pj 

re  = radius  of  the  Bjerhammar  sphere 

rj,2  = rj2  + r 82  - 2 r}  rB  t 

t = cos  il)  , 4>  = geocentric  angle  between  P t and  Pj 

d S = surface  element  at  the  internal  sphere 

S = 4 tt  re3 


Note.  The  gravity  anomaly  (Ag3)  at  the  point  P3  is  defined  as  the  difference 
between  the  real  gravity  at  Pj  and  the  theoretical  gravity  at  a point  p3',  which 
is  located  along  reference  ellipsoid  normal  through  P3  with  theoretical  potential 
at  Pj'  equal  to  the  geopotential  at  Pj.  (See  Heiskanen  and  Moritz,  1967,  p.  83.) 

Due  to  the  irregular  topography  of  the  earth  there  is  rigorously  no 
solution  of  (1.1).  However,  for  the  Bjerhammar's  problem  (discrete  boundary 
value  problem)  it  is  always  possible  to  find  a fictitious  function  Ag*  that  fits 
(1.1)  for  the  observations.  Once  Ag*  is  determined  it  may  be  applied  in 
Poisson's,  Stokes'  or  Vening  Meinesz'  integral  formula  for  the  prediction  of 
gravity  anomalies,  geoidal  undulations  or  vertical  deflections,  respectively. 

It  should  be  stated,  that  in  order  to  be  consistent  with  Stokes'  formula,  the 
harmonics  of  degree  0 and  1 should  be  removed  from  (1.1). 

The  method  of  least  squares  collocation  was  introduced  in  the  principal 
presentations  by  T.  Krarup  (1969)  and  H.  Moritz  (1970,  1972).  In  general  this 
method  is  defined  as  consisting  of  two  parts:  a least  squares  adjustment  and  a 
Wiener  prediction.  In  this  study  we  will  consider  only  the  latter  part  of 
collocation.  If  the  relevant  covariance  functions  are  known,  any  physical 
quantity  (v  3 ) in  an  arbitrary  position  (P3)  on  or  outside  the  surface  of  the 
earth  may  be  estimated  in  an  optimal  way  from  a finite  number  of  observations 
(Ag)  by  the  following  formula: 


(1.2) 


v i = c i (C  + D)  1 Ag 


where 


Cl  = COV  (vlt  Ag) 

C = cov  (Ag,  Ag) 

D = error  covariance  matrix. 


In  particular  we  predict  gravity  anomalies  by: 

(1.3)  Ag!  = c j (C  + D)"1  Ag 

where  the  elements  of  c ! and  C are  given  by  the  spatial  covariance  function  for 
the  gravity  anomalies: 

00 

(!•  1)  CtJ=y  Cn  (re2/r  , rj)n+2  Pn  (cos  0,j) 

n =0 


where  cn  are  the  so-called  degree  variances  of  the  gravity  anomalies  defined 
according  to  the  definition  in  Heiskanen  and  Moritz  (1967,  p.  259): 


(l*5)  cn~  37 r jJ  Agif  d CT 

where  Agn  is  the  anomaly  Laplace  harmonic  on  the  Bjerhammar  sphere. 

The  purpose  of  this  paper  is  to  study  the  relations  between  some  of 
Bjerhammar's  solutions  and  collocation.  We  start  with  solving  Poisson's 
integral  equation  (1.1)  in  different  ways  according  to  Bjerhammar. 


2.  Bjerhammar  Solutions  at  the  Internal  Sphere 

In  the  numerical  application  of  (1.1)  Bjerhammar  used  a finite  set  of 
blocks  on  the  internal  sphere  with  a constant  value  of  Ag*  over  each  block.  In  this 
way  the  following  matrix  equation  is  obtained  from  (1. 1)  for  m observations  and 
N surface  blocks; 


(2.1)  A Ag*  = Ag 

■ N N 1 ml 


-3- 


where 


(2„  2)  Afa  - 


vector  of  observed  gravity  anomalies 


vector  of  unknowns 


coefficient  matrix  with  elements 

_ 3 2 

Ti _rB  r r ds 


AS  k = the  surface  of  block  k 


Applying  the  mean  value  theorem  of  integral  calcuk  s formula  (2.2)  may 
be  rewritten: 


(2.2a) 


2 3 

rl  ~ r« 
4 tt  r„ 


where  k refers  to  a certain  mean  value  point  Pk  inside  ASk.  In  practice  we 
may  approximate  Pk  by  the  center  of  the  block  k.  If  the  matrix  A has  full 
rank  there  is  always  a solution  of  (2. 1)  for  N a m.  For  N = m we  obtain  the 
unique  solution: 


Ag*  = A 1 Ag 


where  A 1 is  the  Cayley  inverse  of  A. 

For  N > m the  solution  is  not  unique  unless  an  additional  condition  is 
satisfied  (condition  adjustment.  Section  2.2).  Finally,  for  N < m and  full  rank 
of  A the  least  squares  solution  minimizes  the  square  sum  of  the  residuals  (see 
Section  2.  3). 


2.1  Generalization 

By  expanding  the  Poisson  kernel  into  a series  of  Legendre's  polynomials, 
formula  (1. 1)  becomes; 


ICC  o / Tb  \n+  2 

= j JJ  Ag*  ^ (2n+l)  (—  J Pn  (cos  i/j)  d S 


-4- 


We  may  also  expand  Ag*  into  a series  of  Laplace's  harmonics  (we  assume  that 
the  expansion  is  convergent): 

CO 

Ag*  = ^ Agn 

n =0 

Due  to  the  orthogonality  of  the  spherical  harmonics  when  integrated  over  a sphere: 

JJ  Agn  P„'  d S = 0 for  n^n 

(2.4)  may  be  rewritten  in  the  following  way: 

(2.5)  Agj  = JJU*^  /(2n+l)  c * ^'j  pn  (C0s  |/>)  d S 

n = 0 


(2.6) 


Agn 

Sc*/(2n+l) 


In  contrast  to  the  degree  variances  in  collocation  (cn),  which  are  a priori  defined 
by  formula  (1.  5),  the  parameters  c n*  introduced  in  formulae  (2.  5)  and  (2.  6)  of  the 
generalized  Bjerhammar  method,  are  more  or  less  arbitrary.  Formally  the 
conversion  from  (2.4)  to  (2.5)  is  valid  for  any  c * such  that  (2.6)  and  the  kernel 
of  (2.  5)  converge.  However,  from  a practical  point  of  view  it  will  hardly  be 
advisable  to  solve  (2.5)  for  a u*  that  is  less  smooth(with  less  attenuating  higher 
degree  terms)  than  Ag*.  Already  the  solution  for  Ag*  is  questionable  in  the 
continuous  case.  See  Section  5,  Molodensky  et  al.  (1962)  and  Pick  (1965).  In 
the  discrete  case  a solution  is  always  possible.  By  changing  c„*  a variety  of 
solutions  to  Bjerhammar's  problem  are  obtained.  These  solutions  for  u*  are 
completely  in  accordance  with  the  solutions  for  Ag*  in  (2.1)  - (2.2). 

Example.-  The  inverse  Stokes  function.  We  insert: 


/c~*  = / 2n  + l 

rB 


-5- 


into  (2.5)  and  (2.6)  (with  lower  limit  n = 2 of  the  summation).  Then  we  obtain 
(for  d S = rB2  d O): 


where 


and 


Ag3  = JJ  M (rJf  il>)  u*  da 

l r / rB\n+1 

M <rJt*)  = l <2n+1>  (“-1)  ("F7 J Pn  (cos  ^ 

3 n=  a 


n = a n =a 


where  we  have  introduced  the  notation  (Heiskanen  and  Moritz,  ibid.,  p.  97): 


Tn  = rB  Aga/(n-l) 


Tn  being  the  Laplace  harmonic  of  degree  n of  the  (fictitious)  disturbing  potential 
T*  at  the  Bjerhammar  sphere. 

M may  be  written  in  a closed  form  (Sjoberg,  1975,  p.  107): 


where 


M ( r j , 0 ) 


s r 3(1  -sV  t ss  - 5 

8rrrj  L jt5  f3 


+ 2 


-i 


j l2  = 1 - 2 s t + s2 

s = rB/rj 

t = COS  <l> 


Note.  In  Sjoberg  (ibid.)  the  term  n = 0 is  included  in  M. 
Further  examples  of  u*  are  given  in  Table  2.1. 


2.2  Minimum  Norm  Solutions 


As  previously  stated,  a matrix  equation 
(2.7)  A u*  = Ag  , m < N 

a N N 1 a 1 

does  not  have  a unique  solution  unless  an  additional  condition  is  imposed.  We 
may  define  this  solution  as  the  one  satisfying  the  following  condition: 


N 

(2.8a)  Uu*|l2=  ~ Y (Uk*)aAS|,  = minimum 

O — J 

k=  1 


or,  with  matrix  notations 

(2.  8b)  (u*)T  Q-1  u*  = minimum 

where 

AS  ^ 

Q"1  = - ( ASa 

s 1 - ASn 

From  least  squares  adjustment  we  obtain  the  following  solution  of  (2.  7)  with  the 
condition  (2.8)  (see  for  instance,  Bjerhammar,  1973,  Ch.  12): 


(2.9)  u*  = Q At  (AQAt)_1  Ag 

with  the  minimum  norm 

(2.9a)  (u*)T  Q_1  u*  = AgT  (AQA1)'1  Ag 

Subsequently,  different  solutions  are  obtained  by  changing  the  number  of 
blocks  at  the  internal  sphere.  In  the  original  paper  by  Bjerhammar  (1964)  most 
studies  were  restricted  to  the  non-singular  case  (N  = m).  Condition  adjustment 
(N  > m)  was  treated  in  1968  and  1969.  It  is  obvious  that  these  solutions  from  (2.9) 
are  generally  more  cumbersome  to  compute  than  (2.3).  However,  the  computa- 
tional effort  is  drastically  reduced  to  that  in  the  non-singular  case,  if  we  let  N 
approach  infinity  for  well-behaving  surface  elements.  This  we  show  next. 

-7- 


' 


We  rewrite  (2.9)  in  the  following  way: 


(2.10) 

where 

(2.11) 


u*  = QAt  X 

X * (AQAt)_1  Ag 

N 

(AQA’ju  = S £ (A)ik  (A)Jk  ASk"1 


The  coefficients  (A)Jk  are  given  by  the  generalization  of  (2.1)  and  (2.2a)  by 
(2.  5): 


(2.12) 


(A)jk  = S I /<2n+1)  c»*  (““)  P*  (cos  0)  ASk 

n=0 


By  letting  N go  to  infinity  in  such  a way  that  the  largest  diameter  of  all  ASk 
approaches  zero  it  is  shown  in  the  Appendix,  Proposition  A.l,  that  (AQAr),i 
becomes.* 


(2. 13) 


r"*  / v \n*  ^ 

Cm  = lim  (AQ  Ar)u  = 2,  ct  {^T~p)  Pn  (cos  0,,) 

N — oo  a = 0 r*J 


Thus  we  obtain  in  the  limit; 


(2. 13a)  X * C-1  Ag 

where  the  elements  of  C are  given  by  (2. 13).  Applying  (2.  7)  for  the  prediction  of 
new  anomalies  from  the  original  set  Ag  we  obtain: 

Ag, = c, X 
or 

(2.14)  Ag,  = c,  CT1  Ag 


-8- 


where 


c i = lim  A 1 Q AT 

N “•  00 

The  elements  of  ct  are  also  given  by  (2.13).  Furthermore  it  follows  from  (2.9a) 
that  the  minimum  norm  in  the  limit  is  given  by: 

(2.15)  ||u*||3  = AgT  C"1  Ag 

We  notice  that  (2. 14)  and  (2.13)  are  exactly  the  solutions  obtained  with  the 
collocation  formulae  (1.3)  -(1.4)  with  cn  substituted  by  cn*  . Hence,  for  each 
type  of  degree  variances  (cn*)  in  collocation,  there  is  an  identical  minimum  norm 
solution  in  the  generalized  Bjerhammar  approach.  Some  of  these  relations  are 
given  in  Table  2. 1.  The  derivations  are  given  in  Sjoberg  (1975). 


Table  2. 1 

The  Relations  Between  the  Degree  Variances  in  Collocation 
and  Some  Minimum  Norms  in  the  Generalized  Bjerhammar  Method' 


/c„*/(2n+l) 

Norm 

Explanation  of  Symbol 

1 

Ag* 

gravity  anomaly 

(n-l)/rB 

T* 

disturbing  potential 

4 n(n-l)/(2n+l) 

<P* 

density  layer 

4nn(n-l)/rB(2n+l) 

M * 

double  layer 

(n-l)/(n+l) 

6* 

gravity  disturbance 

y(n-l)//n(n+l) 

e* 

deflection  of  the  vertical 

(n-l)//(2n+l)(n+l) 
t*  1 n!  (n-l)/(n+t') ! 

6* 

d(t,)T* 

* rv 

potential  gradient 

the  i/th  derivative  of  T 

y = normal  gravity  at  the  reference  ellipsoid 


t from  Sjoberg,  (1975). 


-9- 


r 


Subsequently  there  is  a duality  between  the  choice  of  degree  variances  in 
collocation  and  the  minimum  norm  in  the  generalized  Bjerhammar  technique. 
Such  an  approaohwith  a minimum  norm  was  also  emphasized  by  Krarup  (1969). 

Remark,  c „ of  formula  (1.4)  was  replaced  by  (2n+l)  c„  in  Sjoberg  (ibid.)  and 
by  (2n+l)  a„a  by  Lauritzen  (1973,  p.  73)  and  Bjerhammar  (1974 and  1977a, b). 

See  also  Krarup  (ibid.). 

Also  in  the  generalized  Bjerhammar  technique  the  observation  errors 
can  be  taken  into  account  (cf.  formula(1.2)J.  In  this  case  we  start  with  the  model: 


where 


A u*  = Ag  - € 


E { e cT)  = D 


Rearranging  this  equation  as: 


...f] 


= Ag 


we  readily  obtain  the  solution: 


[AQAr  + Df1  Ag 


which  minimizes 


(u*)T  Q'1  u*  + CT  D"1  ( 


The  predictions  are  given  by 

(2. 16)  v,  = A , Q At  (AQ  AT  + Df1  Ag 


and  for  N-®  this  solution  and  (1.2)  are  identical.  In  this  case  the  minimum 
norm  becomes: 


(2.17)  Agf  C'1  Ag  + fT  D"1  c 

-10- 


(2.9b) 


1 r rT 
u*  Q A 

( ~ D 


2.3  Least  Squares  Solutions  for  u* 

If  we  apply  formula  (2.7)  with  m > N,  this  matrix  equation  will  not  in 
general  be  consistent.  Adding  a residual  vector,  we  obtain  the  following  model: 

(2. 18)  A u*  = ig  - ( , m > N , E {c  eT  } = p o2 

a N N 1 a 1 a 1 

This  system  may  be  solved  with  ordinary  least  squares  adjustment  by  elements 
with  the  solution: 

(2.19)  u*  = (A'PA)'1  At  P Ag 

which  minimizes  the  square  sum  of  the  residuals  (fT  P ().  The  variance  of  unit 
weight  (a2)  is  estimated  by: 

(2.  20)  s2  = AgT  P (Ag  - A u*)  /(m  - N) 

and  the  covariance  matrix  of  u*  after  adjustment  is  estimated  by: 

(2.21)  Quu  = s2  (AT  P A)"1 

All  these  derivations  follow  from  elementary  least  squares  adjustment  (see  for 
instance,  Bjerhammar,  1973,  Ch.  11). 

The  solution  u*  can  now  be  used  for  the  prediction  of  derived  quantities 
(v ,)  in  an  arbitrary  position  P,: 

(2.22)  v,  = A t u* 

where  A t is  the  operator  (vector)  that  relates  v,  to  u*. 

The  prediction  variance  may  be  determined  in  the  following  way: 

e,  = v,  - v,  = v,  - A,  (u  + cu) 
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€ I = vf  + A , (uuT  + eu  cuT)  A,'-  2 A ! (u  + tu)  v,T 


Taking  the  expectation  of  both  members  we  obtain  the  prediction  variance: 
m i = (7y  i + A j (C  (in  ^ Q uu  ^ ) A i - 2 A i C vjv 

* 1 

where 

< = E{vt2}  | 

C uu  = E { uuT  } ! 

r r, 

Cuv  = E {uvM 

Furthermore  we  have  assumed  that  E [( u viT ) = 0.  The  prediction  variance  may 
be  written: 

(2.23)  m 2 = 0s  A t Q uu  A J + m/ 

where 

m*  — (T  »i  + AjCuu  A|  - 2 AjCu* 

The  first  term  of  (2.  23),  which  is  caused  by  the  error  of  u*,  is  easily  determined 
in  the  adjustment.  The  second  term  (mA2)  is  due  to  the  error  of  A , . The 
determination  of  m 2 requires  that  the  covariance  relations  are  known. 

The  least  squares  solutions  are  favorable  first  of  all  when  a large  number 
of  observations  (m)  are  available.  It  should  be  noted  that  the  solutions  for  m > N 
do  not  include  pure  prediction.  A filtering  of  the  data  is  always  present  in  this 
type  of  prediction. 

2. 4 Impulse  Approach 


Bjerhammar  (1974)  introduced  a method,  where  the  unknowns  (u*)  are 
located  in  discrete  points  at  the  internal  sphere.  This  method,  called  the  Dirac 
method,  gives  the  following  model: 


(2.24) 


where 


u*  (r)  = £ u\  6 (r-  rt 


u\  = 


(r & , <p,  A.)  = coordinates  of  a current  point  at  the  internal  sphere 

(rb,0k,Ak)  = coordinates  of  a selected  carrier  point  Pk  at  the 

internal  sphere 

unknown  associated  with  Pk 


N = number  of  unknowns 


6( ) = Dirac's  delta  function,  defined  by 


(2.24a) 


Jju*  (r)  6 (r  - rk)  d o = u*  (rk) 


Inserting  (2.24)  into  (2.5)  we  arrive  at 


(2.25) 


where 


(2.26) 


ni 

Agj  = £ uk*  K (?j , rk) 


CD 

^ r0  n + a 

K (rJt  rk)  = y /(2n+l)  c „*  Pn  (cos 


For  m observations  (Agj)  (2.25)  gives  an  exact  matrix  equation  which  may  be 
solved  by  means  of  condition  adjustment,  least  squares  adjustment  or  direct 
solution  dependent  on  whether  m < N,  m > N,  or  m = N.  A generalization  of 
the  method  is  obtained  if  we  also  allow  for  carrier  points  located  outside  the 
Bjerhammar  sphere  (reflexive  prediction,  Bjerhammar,  1974). 

It  is  obvious  that  this  method  is  a generalization  of  the  commonly  used 
buried  mass  point  method.  In  fact,  the  latter  is  obtained  in  the  special  case 
(cf.  Table  2. 1). 

ST*  = const.  (n-l)//(2n+l) 


i * * j"-  ->  ^ •• 


The  non-singular  Dirac  method  with  carrier  points  located  at  the 
intersections  of  the  internal  sphere  with  the  radius  vectors  of  the  observations 
is  very  similar  to  collocation.  The  auto-covariance  function  for  Ag  in  colloca- 
tion (with  internal  sphere  radius  re): 


(2.27)  C(j,i)=£  cB(rB3/r1r1)a+2P»(cos^n) 

n = a 


corresponds  to  the  following  kernel  function  in  the  Dirac  method: 

00 

(2.28)  C (j,i)  = £ y(2n+l)  c * (r0/r1)n'''2  Pn  (cos  ib^) 


n=3 


where  r0  is  the  selected  internal  sphere  radius  (not  necessarily  the  same  as  rB). 
Hence,  C (j,  i)  of  (2.27)  and  (2.28)  are  identical  for: 


and 
(2.  29) 


c>  = c * = 2n+l 


r„  = rBs/ri 


Consider  the  case  rt  = r = radius  of  the  mean  earth  sphere: 

r0  = r - h0 


and 


rB  = r - hB 


where  h0  and  hB  are  the  depths  of  the  internal  spheres  from  the  mean  earth 
sphere.  By  inserting  these  expressions  for  r0  and  r&  into  (2.29)  we  obtain: 


ho  = 2 hB  - hBa/r 


or  (after  omitting  the  last  term) 

(2.30)  h o = 2 hB 
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1 


Thus  we  have  shown  that  in  the  special  case  c n = c * = 2n+l  and  the  radii  of 
the  observation  points  are  constant  (r),  the  kernel  function  C (j,  i)  of  the  Dirac 
method  and  collocation  are  identical  whenever  the  depth  to  the  Bjerhammar 
sphere  in  the  previous  method  (h0)  is  twice  that  in  the  latter  method  (hB).  This 
result  holds  also  for  the  predictions  of  the  two  methods  (see  Secfions  5 and  8). 

See  also  Bjerhammar  (1977,  a,b). 

From  numerical  point  of  view  it  is  of  importance  that  the  kernel  function 
(2.28)  is  unsymmetric  in  contrast  to  (2.27). 


3.  Stability  of  the  Solutions 

It  is  of  great  importance  for  the  prediction  results  that  the  matrix 
equations  are  well  conditioned.  In  this  section  we  are  going  to  compare  the 
stability  of  the  following  systems  in  different  cases: 


(3. 1) 

A u* 

= Ag 

a a 

a 1 

and 

(3.  2) 

c X 

= Ag 

a a 

a 1 

These  equations  were  introduced  in  (2.7)  and  (2. 13a).  Formulae  (3. 1)  (m  = N) 
and  (3. 2)  (N  = ® ) may  be  regarded  as  the  extreme  cases  of  (2.  7).  From  the 
derivations  in  Section  2.2  we  may  also  consider  (3.2)  as  an  intermediate  system 
of  equations  in  collocation. 

A suitable  measure  of  the  stability  of  a matrix  system  is  the  condition 
number  defined  as: 

(3»  3)  x = A a,*  / A B ln 

where  A refers  to  the  eigen  values  of  the  coefficient  matrix  of  the  system  (A  and 
C).  The  condition  number  has  the  following  bounds: 


1 s x £ 00 

where  the  lower  bound  is  the  ideal  situation  and  00  is  a completely  unstable 
(singular)  system.  Our  main  problem  is  therefore  to  determine  the  eigen  values 
A ■ in  and  A a*x  • 
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We  start  with  a two-dimensional  example.  We  assume  that  Ag  is 
observed  at  m regularly  distributed  points  on  a circle  of  radius  r.  The 
unknowns  (u*)  are  located  on  a circle  of  radius  r&  . Applying  the  impulse 
approach  [formula  (2.25)  with  m = N]  the  system  of  equations  becomes: 

N 

(3.4)  Y K (6j-  0k)  u„*  = Agj  ; j = 1,  2 N 


where  K (0j  - 0k)  is  the  two-dimensional  kernel  function  corresponding  to  u*  and 


0 p = p 2 rr  /N 


The  Fourier  series  of  K is: 


where 


K (0)  = Y bB  eln0 


b„  = 2^  J K (6)  e"ln6d  6 


Furthermore,  the  eigen  values  (A 4)  of  (3.4)  are  given  by: 


(3.7)  V K (0j-0k)  x[£)=  A£  x}n  ; i,  j = 1,  2 N 


where  X/  ' is  the  j th  element  of  the  £th  eigen  vector.  Let  us  try: 


X ]'^=  cos  x = (eu+  e lx)/2 


where 


Using  the  relation 


x = j i 2 tt/N 


Ye‘^<*-q)  = { 


N if  (,=  q + Np  ; p = 0,  ± 1,  ±2,  . . . 
0 otherwise 


we  can  easily  arrive  at: 


l Xj(V(q)=  {N/2  if  *=q 

^ ^ 0 otherwise 


Thus  (3.  8)  satisfies  the  necessary  property  of  eigen  vectors  being  orthogonal  to 
each  other.  We  have  also  shown  that  (3.  8)  should  be  normed  by  /2/ N in  order  to 
be  the  eigen  vectors  of  (3.  7).  Finally,  X must  satisfy  (3.  7)  for  some  eigen 
values  A*.  Inserting  (3.  8)  and  (3. 5)  into  (3.  7)  we  obtain  from  the  left  member: 


S an  JL  //t  .2 17  . 2Tt  « 

£ bn  etnJ  n Vel(^'n)  n * = NeltJ»  J be+NP 


k = l 


Here  we  have  used  (3.9).  Thus  we  obtain: 


= N 


+ NP 


P = —CD 


where  bn  is  given  by  (3.6).  If  K is  Poisson's  kernel,  then  (Seeley,  1966,  p.  14) 


b„  = s,n|  = (rB/r),n| 


and 


>*l 


The  sum  is  given  in  a closed  form  in  the  Appendix,  Corollary  A. 2.  We  obtain: 


H = N (sC  + sN“V(l  " sN)  , 1 = 1,  2,  ....  N 


The  condition  number  is  finally  given  by: 

*i  = (l  + sN)/(s^+  sN-^) 
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where 


0 = j = integer  part  of  N/2 


An  analogous  derivation  of  the  condition  number  of  (3.2)  gives: 

xa  = (1  + S2N)/(S28  + s2N-sS) 

so  that  for  large  N 

Xa/Xx  % s~n/s 

It  is  obvious  that  the  system  (3. 1)  is  more  stable  than  (3.2)  for  a given  depth  to 
the  Bjerhammar  circle. 

Next  we  proceed  to  the  planar  cases  of  (3. 1)  and  (3.2).  We  are  going  to 
use  an  approximate  similarity  transformation  of  the  coefficient  matrix  to  diagonal 
form. 

The  operator  C is  said  to  be  unitarily  equivalent  to  S,  if: 


S = U C U' 


where  Si  j = Xu  6jj  , 6U  is  the  Kronecker's  symbol  and  U and  U-1  are  unitary 
transformations,  inverse  to  one  another.  For  a matrix  it  is  a similarity 
transformation  to  diagonal  form.  According  to  Moritz  (1966)  and  Schwarz  (1971) 
such  a unitary  equivalence  exists  between  c(x,y)  and  its  Fourier  transform: 


(3. 10) 


00  CO 

S (u,v)  = U • c (x,  y)  =*  JJ  c (x,y)  e-  1 ^xu+  yv)  d xd  y 


This  formula  may  therefore  be  used  for  an  approximate  determination  of  the 
eigen  values  of  (3.1)  - (3. 2).  If  m0  and  no  are  the  numbers  of  blocks  in  the 
x-  and  y-  directions,  we  obtain  (cf.  Schwarz,  ibid.): 


(3.11) 


A.  nl  = a s(J^SL  , 

\2n0  + l 2m  0 + 1 / 
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where 


0 ^ n s n0 


0 s m £ m0 


and 

a=  ||c|| 


The  covariance  functions  will  be  used  in  the  planar  approximation.  Furthermore, 
equal  area  blocks  are  assumed  with  n0  = m0. 

We  are  now  going  to  determine  the  condition  number  for  the  Poisson 
kernel.  This  kernel  corresponds  to  cn  = 2n  + l and  the  minimum  norm  of  Ag* 
(see  Table  2. 1).  The  planar  equivalence  to  the  operator  A of  (2. 1)  is: 


(3. 12) 


A = 


JJ  a(x,y) 


dx  dy 


where 


(3. 12') 


a <x.y)  = ~ / (x2  + y*  + h2) 


h = z + b/2 

b = 2 (R  - rB) 

R = mean  earth  radius 

r = radius  of  the  computation  point 

rB  = radius  of  the  Bjerhammar  sphere 

z = r - R 

Formula  (3. 12  ) has  the  following  Fourier  transform-. 


S(u.v,=  JJ  a (x,  y)  e x P { - i (u  x + v y) } d x d y 


(3.13) 


= 2n  h 


r-  ■■■;  J 0 (oer ) d r = 2nexP  {-hte) 

J / —2  ^ ^2  v*V  2 


•J  / _,2  , i 2.3/ 2 

o (r  +h  ) 
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where 


and 


Hence 


and 


(3.14) 


r = v^x2  + y3 


J o (to  r)  = zero  order  Bessel  function. 


.'=  ^ua+  v3 


An,  = 2tt  a exP  (-2nh>/nli  + m8  /(2n0  + 1)  } 


x«  = ex  P (2n«/2hno/(2n0+l)]*  e..  P {n  '^2  h ] 


The  planar  approximation  to  the  Poisson  covariance  function  is 
(Moritz,  1976,  p.  41): 


(3. 15) 


C (x,y)  = C (r)  = 


B/Z3 


(l  + rVz3) 


,a  /™a.3/a 


where  Z = z , + Zj  + b and  B = 2 Ra. 

The  Fourier  transform  of  the  covariance  function  is : 
(3. 16)  S (uj)  = 2 tt  B ex  P [-Z  u) } 

Thus  we  obtain: 


xc  = ex  P [2rr^2  Z n0/(2n0+  1)} 


Assuming  that  z , - z ) ^ z we  have  Z =»  2h  and 


(3.17) 
so  that 

(3. 18) 


Kc=exP  {4/TrTh  n0/(2n0+  1)}  as  ex  P [2v/2nh  } 


xc  ^ x. 
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From  formula  (3. 18)  it  is  obvious  that  the  condition  number  (x»)  of  the  "original" 
equation  (3. 1)  is  squared  by  letting  the  number  of  blocks  on  the  internal  sphere 
approach  infinity  (N  ) (and  then  solving  the  system  (3.  2)).  A considerable  loss 
of  stability  is  achieved.  Again  we  emphasize  that  the  equation  system  (3. 1)  is 
generally  unsymmetric  in  contrast  to  (3.2),  which  fact  has  to  be  considered  in 
the  numerical  solution  of  the  system.  It  is  not  recommended  to  form  normal 
equations  in  (3. 1)  in  order  to  obtain  a symmetric  system,  because  this  procedure 
will  increase  the  condition  number  to  that  of  system  (3.2).  (3. 1)  should  be 
solved  directly,  for  instance,  by  Gauss  elimination  or  in  an  iterative  way  by 
successive  approximations. 


3. 1 The  Effect  of  Smoothing  (Noise) 

We  assume  that  the  matrices  A and  C of  formulae  (3.1-2)  are  substituted 
by*. 


A = A + D 


and 


C = C + B D 


where  the  elements  of  D are  given  by*. 

k if  i = j , k > 0 

du  = { 

k 0 otherwise 


and  B = 2 R°  {see  formula  (3. 15)  1 . 

Then  A and  C correspond  to  the  kernel  functions: 

a (x,y)  = a (x,y)  + k 6 (x,y) 

and 


with  the  spectra 


c (x,y)  = c (x,y)  + B k 6 (x,y) 


S.  (u,  v)  = S.  (u,v)  + k 
Sc  (u,  v)  = Sc(u,v)  + B k 
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and  the  eigen  values 


Xw  . 
A na 

>) 

X n« 

a k 

-(c) 

x„.  = 

(c) 

XB«  + 

a B k 

Here  6 is  the  delta  function  and  S,  A.  and  a are  the  same  as  in  the  previous  section. 
In  this  case  we  obtain  for  the  Poisson  kernel  [cf.  (3. 14)  and  (3. 17) } : 


and 


X,~  k +ii  exp  {-ttIwS  } <exp{/*TTh} 


k + 2 TT 

k + 2nexpf-2y5h}  <expf2/2rrh} 


We  notice  that  the  smoothing  stabilizes  both  systems  of  equations  (cf. 
Section  6.3). 


4.  On  the  Convergence  of  an  Iterative  Solution 

In  this  section  wc  arc  going  to  investigate  a condition  for  convergence  of 
an  iterative  solution  of  the  matrix  equation  (3.2): 

(4. 1)  C X = Ag 

Using  the  ordinary  method  of  successive  approximation,  we  obtain: 

(4.2)  X(k)=  Ag  + (I-C)  X(k_l);  k = 1,  2,  ... 
and 

AX(k)=  X(k)-  X(lt_l)=  (I-C)(X(k_l)- X(k_8))  . 

(I-C)k_1  (X(l)-X(0)) 
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Hence 


|AX(y)||  s ||l- C If"1  ||X(1)-X(0)| 


A sufficient  condition  for  the  convergence  of  (4.2),  i.e. 


is  therefore 


lim  J AX  | = 0 


I - C I < 1 


max  Y 1 6 u - c u | < 1 


I 6u-  c,j  | = { 


|l  - c u | , if  i = j 

|c  u I • if  i t j 


so  that  the  condition  may  be  written: 


(4.3a) 


(4.3b) 


a 

Y |c,j  | < 2 , if  clt  > 1 


n > Y |c  u | , if  c u 


For  finite  Cn  we  may  always  divide  each  row  of  the  equation  system 
with  a sufficiently  large  number  to  ensure  that  the  diagonal  elements  of  C are  equal 
to  or  less  than  unity.  Hence  we  can  limit  our  study  to  the  condition  (4.  3b). 

First,  we  study  this  condition  for  Poisson's  integral  equation  for  the 
circle  (Sealey,  1966,  p.  14): 


• , * » *•  i _ •*%  *•  9 


(4<4)  Ag  j = J*  a (j,  i)  Ag*  d 0, 

- TT 

where 


a (j.‘) 


1 1 - s2 

2n  l+s2-  2s  cos  4>  u 


1_  1 - sa 

2n  |s-e^,J|2 


s = rB/rj  (<  1) 

= 8j-  0, 


Suppose  that  we  apply  (4.4)  in  a discrete  approach  with  equal  spacing  between  the 
N unknowns  (Ag**).  We  obtain: 


(4*  5)  A Ag*  = Ag 

1 N N 1 ■ 1 

where 

(4.6)  A Jk  = a (j,k)  A6  = (1  - sa)/N  |s  - e^®3  ~ |s 
A 6 = 2 tt/N 

0k  = 2 nk/N  ; k = 1,  2,  . ...  N 


In  the  non-singular  Dirac  approach  (N  = m)  we  arrive  at  the  following  condition  from 
(4.3b)  and  (4.6)  (for  0J=  0): 


or 

(4.7) 


_2_  l-s*  > J_  f 1-  sa 
m (l-s)a  m L |s_e i0 k |a 


S (m,s)  = — t — — — iUlSL  < o 

m - |s_e,e’'|2  m(l-s) 

k=  l ^ * 


The  Poisson  kernel  has  the  following  Fourier  series; 
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The  corresponding  relation  for  collocation  is  obtained  by  substituting  s by  s2  in 
the  previous  formulae.  The  results  are  shown  in  Table  4.1. 

For  comparison  we  also  derive  the  relation  between  m and  s in  the  case 
of  condition  adjustment  (N  > m).  Then  we  arrive  at  the  following  condition  from 
(4.3  b): 


S (m,  N,  s)  = max  Y (AQAr),j  -2(AQAt),i  < 0 

i 

j=i 

where 

N 

(A  Q A )ij  = N V A ik  A jk 

k=i 
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Table  4. 1 


For  large  m we  finally  arrive  at: 


m < m o = 2 


These  are  the  limiting  relations  between  s and  m when  using  collocation  [cf. 

(4.8  a-b)]. 

What  Is  the  relation  between  s and  m that  satisfies  the  condition  (4.3  b) 
for  the  generalized  integral  equation: 


(4.4) 


where 


4gj  “ J a (j,k)  u*  d 0* 


s"V(9>-e»> 


Following  the  previous  derivations,  we  obtain  in  the  non-singular  case  of  the  Dirac 
method  (cf.  (4.7)]: 


00  CO 

S (m,  8)«£  /5SS  SIJI  * - ^ j s'n*  < 0 


Thus  the  sufficient  relation  between  s and  m for  the  convergence  according  to 
(4. 3b)  is  very  much  dependent  on  c*  . 

Example;  /c*  = |n|  + l 

S(m,s)=|  |j|m  sIJ"  + £ sIJ"-  2 (£  |n|s""-£  s'n|)  = 


\ ds  / \ j_s»  m m 1-8  / 


- 2m  -i— - ± _S_t  l+  -i.i.  -£_<  o 

(1-s’)3  m (l-s)a  1-8  m m 1-8 


For  large  m we  have  approximately 


. 4s  4s  _ . 
m - 2 - < 0 


(1-s)3  1-s 


m < m o = 2 + 


s > So  = 1 


4s (2-s) 


2 


For  m0-®,  s0-  1 and  h/a-  2/tt  . Again  we  obtain  2n  as  a sufficient 
ratio  between  the  height  of  the  observations  and  the  spacing  of  the  observations 
for  convergence  in  the  continuous  case  (cf.  Table  4.1). 

We  now  investigate  Poisson's  integral  equation  in  the  planar  approximation 
(see  formula  (3.12)).  We  assume  that  a square  of  side  B is  divided  into  N surface 
elements  of  side  b = B/N  , where  N =</n" . Then  we  obtain  the  following  coefficients; 


A - 


2n[(xrxk)3+(yJ-y1,)3+h3]3/8 


= c/2tt  ((j)t  - kjt)3+  (jy  -ky  )3  + ca)3/3 


where 
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In  the  non-singular  case  (N  = m')  the  condition  (4. 3b)becomes  (for  j*=  jy  = m'/2): 


S (m,c)  = 7 t — — — — —/s  ~ ~ < 0 

2 - - [(m  /2  -k  x)2+(m  /2  -kr  )2+  c2]3'2  nc3 

kx»iky=l 


The  critical  relation  between  m'  and  c is  shown  in  Table  4.2. 


Table  4.2 


The  Ratio  c=  h/b  Satisfying  S(m,c)  = 0 as  Given 
for  Different  m.  b=  block  size,  h=  height. 


nr/ 

Dirac  Method 

Collocation 

50 

0.57 

0.285 

100 

0.54 

0.  270 

500 

0.52 

0.260 

1000 

0.52 

0.260 

5000 

0.51 

0.255 

10000 

0.  50 

0.250 

The  result  of  the  computations  is  that  the  ratio  h/b  should  be  less  than 
0.5  for  an  uncritical  convergence  of  the  non-singular  Dirac  solution  with  Poisson's 
kernel  (see  also  Koch,  1968).  If  the  method  of  collocation  is  used,  the  ratio  should 
be  less  than  0.  25.  In  this  study  we  have  based  the  results  on  the  condition  (4.  3b) 
fora  solution  by  successive  approximations.  However,  Schwarz  (1971)  emphasized 
that  more  generous  ratios  (h/b)  may  be  allowed  by  using  other  numerical  methods. 
The  figures  given  in  this  section  should  therefore,  first  of  all,  be  used  to  compare 
the  stability  of  different  methods  (Dirac  method,  condition  adjustment,  collocation) 
with  each  other. 

In  conclusion,  the  Dirac  method  with  (m=  N)  gives  a more  stable  solution 
than  collocation  for  a given  radius  re . In  the  case  of  condition  adjustment  (N  > m) 
the  stability  of  the  Dirac  method  is  roughly  the  same  as  for  collocation.  This  study 
was  restricted  to  the  Poisson  kernel  functions,  and  the  conditioning  changes  with 
the  choice  of  degree  variances  (c  „ and  cn*).  In  general,  however,  the  above 
tendency  can  be  expected  for  an  arbitrary  type  of  kernel  function. 
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5.  Predictions  According  to  Bjerhammar 


Using  the  collocation  formula  (1.2)  any  geophysical  quantity  vt  can  be 
predicted  from  a vector  Ag  of  gravity  anomalies  by: 


(5.  la)  v i = c i X 

where 

(5.1b)  X = (C  + D)-1  Ag 


In  the  generalized  Bjerhammar  methods  the  unknowns  X of  collocation  are  replaced 
by  the  unknowns  u*  at  the  internal  sphere  (see  section  2).  The  prediction  of  v t 
becomes  accordingly: 

(5.  2)  v i = Fi  u* 


or 

(5.  2) 


N 

v i = y Fiu  u£ 

k=j 


where  F ( is  the  relevant  coefficient  matrix  (with  elements  Fix)  relating  u*  to  vt  . 
Formula  (5.2)  is  the  discrete  form  of  the  integral  formula; 


(5.2") 


r‘=  T JJf“u*ds 


and  the  kernel  function  f u,  and  the  coefficient  Fi*  are  related  by  {cf.  formulae  (2.2) 
and  (2.  2a)  J : 

(5.3)  Fix  = Jj  fix  d S = f,x  ASk/S 

ASx 

where  k corresponds  to  some  mean  value  point  inside  ASx.  In  practice  this  point 
is  approximated  by  the  center  of  ASx  • The  function  f tx  can  be  derived  from 
formula  (2.  5)  for  each  specific  quantity  v t.  For  v i = Agt  we  have: 


(5.4)  fix  = Y •/(2n+l)c?  (re/ri)n+3  Pn  (cos^ix) 

• = o 
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and  File  = Aik  according  to  (2. 12).  Further  examples  are  given  below. 


d 


where  y is  the  theoretical  gravity  at  the  reference  ellipsoid.  In  the  special  case 
ri=  re  and  cn*=  2n+l  (u*=Ag+)  (5.7)  inserted  into  (5.2  ) yie  ds  Stokes'  formula 
(Heiskanen  and  Moritz,  ibid.,  p.  94). 


Prediction  of  Deflections  of  the  Vertical 

The  two  components  of  the  deflections  of  the  vertical  (4, T7)  are  related 
to  the  height  anomaly  C according  to  ( Heiskanen  and  Moritz,  ibid.,  p.  312); 


where  <pt  and  are  the  geodetic  latitude  and  longitude.  By  using  the  following 
relations  from  Heiskanen  and  Moritz  (ibid. , p.  113)  and  Abramowitz  and  Stegun 
(1964,  p.  334): 


3 th  d lb 

— — = - cos  a — r-  = - sin  a cos  <p 

d<p  dX 


and 


3Pn(cos  Ip)  _ 
dip 


- sin  ip 


d Pn(cos^) 
d cos  Ip 


n (n+1) 
(2n+l)sin  <l> 


Pn-i;cos 


where  Of  is  the  azimuth  from  the  fixed  (computation)  point  to  the  moving  point,  we 
obtain  the  following  kernel  functions  for  £ j and  r)t  : 

r COS  01  lk  . 

(5.8)  f -jc  = \ . \ x 

tsin  a tK  J ysin^ik 

X V ( — ) {P n + l (COS  <i))  — P n-i  (COS  (P)  ] 

n-1  v 2n+l  \rt/ 

n=2 

where  cos  Of  jj(  refers  to  £ i and  sin  ot  uc  to  Tj  Formula  (5.8)  inserted  into  (5.  2 ) 
gives  a generalization  of  Vening  Meinesz'  formula  (Heiskanen  and  Moritz,  ibid. , 
p.  114). 

Prediction  of  the  Vertical  Gradient  of  Gravity 

The  kernel  function  for  vt=  3Ag[/3ri  is  easily  obtained  from  (5.4).  The 
result  is; 


(5.9) 


00 

fix  = --V  ,/(2n+l)  c f (n+ 2)  (r8/r1)n+2  P„  (cos  l/Jlk) 

ft 


n=  0 

For  certain  sets  of  c*B  the  above  kernel  functions  may  be  given  in  closed  forms 
(cf.  below  and  Sjoberg,  1975,  section  A. 2). 


Note.  In  the  impulse  methods  (section  2.4)  Ft*  = f i* . 


5. 1 Minimum  Norm  Solutions 


The  prediction  of  a quantity  v ( according  to  any  of  the  methods  of 
Bjerhammar  described  in  section  2 is  given  by  formulae  (5.  2)  and  (5.3)  for 
any  finite  number  of  unknowns  (N).  For  example,  in  the  minimum  norm 
solutions  (section  2.2)  the  prediction  formula  becomes  from  (5.2)  and  (2.9): 


vt  = Fi  Q AT(A  Q A7)'1  Ag 

From  this  formula  we  obtain  in  the  limit  (N^^Hcf.  2. 14): 

(5.10)  v , = k , C-1  Ag 

where  the  elements  of  the  matrix  C are  given  by  formula  (2. 13): 


and 


00 

Cij  = Y cn*  ( re/  r i r/  + 2 P„  (cos  Ip  u) 

n=0 


k i = lim  Fi  Q AT 

N -•  oo 


The  following  elements  k11(of  ki  are  obtained  in  accordance  with  formulae  <2. 13), 
(5.4)  and  (5.  6)-(5.9).  The  coefficients  c*  are  given  bv  the  selected  minimum 
norm  (see  Table  2. 1). 


(5.4)  kuc=y  c*„  sJt2pn(t)  ; forvi=Ag( 

n = o 

(5.6)  klk=^Y  ll-sn  + 2Pn(t)  ; for  v i-  Ci 

y -<  n-l 

n =3 


(5.  8) 


k lk 


1 

ysin^Ik 


l 


c * n(n+l) 
n-l 


n=  3 


f cos 

; for  v i = 

{5'} 

^sin 

«lkJ 
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(5.9) 


where 


ku  = - - Y C*a  (n  + 2)  sn+a  P„  (t) 

-i 

n=  0 


; for  v i = 


d 

3 rt 


s = 


tb/t  i 


and  t = cos  4>  ik 


The  above  predictions  from  formula  (5. 10)  are  identical  with  the  predictions  using 
collocation  (without  noise).  Erroneous  observations  can  be  considered  by  adding 
a noise  covariance  matrix  to  C of  (5.10)  before  inversion  {see  section  2,  formula 
(2. 16)}.  The  only  type  of  observations  considered  in  Bjerhammar's  methods  are 
free-air  gravity  anomalies.  However,  as  in  collocation,  heterogeneous  data  might 
be  included  in  the  predictions. 

Finally,  we  mention  that  for  many  norms  the  above  kernel  functions  can 
be  written  in  closed  forms.  As  an  example  we  give  the  kernel  functions  for 
|Ag*,i  = min.  (c  * = 2n+l): 


(5.4  ) 
(5.6") 
(5.8") 


(5.9  ) 
where 

and 


k ix  = s (1  - s )/ 1 ; for  v t = Ag  j 

kix=  rj  s“  {2/£  + 1 -3£-s  t(5+3&i  0/2)};  forv(=Ti 

i s3  . f \ f 2 3 c 3st(l+Z)  . 

y l sin  a,x  J L g3  £ £ 0 

4 


+ 3 in  0/2  } ; for 


k ix=  ; {l  - 5s~  + 3(1  -sa)  /£2};  for  v i=  d Agi/S  r x 

2r  ti 


s = r S' r i r x 


t = cos  lk 


= (1-  2 st+  sa)^  , 0 = 1 — s t + 


More  examples  with  derivations  are  given  in  Sjoberg  (ibid. , Appendices  A.  2 and 
A.  3). 
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5.2  Model  Studies 


In  Sjoberg  (1975,  section  19)  some  minimum  norm  solutions  [formula 
(5.10)}  for  C,  Ag  and  0 (=^2+VS)  are  computed  for  two  simple  earth  models. 
Some  of  these  results  are  reported  below.  For  details  we  refer  to  Sjoberg  (ibid.). 
The  following  norms  are  used; 


Method 

Norm 

c n*/(2n+l) 

1 

Ag* 

1 

2 

T* 

(n-1)5/  re” 

3 

Lauritzen1^ 

(n-l)/((2n+l)(n-2)) 

4 

0 

(n-l)s/(n  (n+1)) 

empirical  norm  from  Lauritzen  (1973,  section  10.3) 

a)(0*)s=  {$*?+  (v*)z 


5.2.1  Bjerhammar's  Model 

Bjerhammar's  model  consists  of  a homogeneous  sphere  with  a spherical 
and  homogeneous  mass  disturbance  M (Figure  2),  where  M = 8.37758  x 10  °kg  and 
a = 6362  km.  The  radius  of  the  main  sphere  is  6370  km. 


Figure  2.  Bjerhammar's  Model.  The  disturbing  sphere  M is  located  with  its 
center  inside  the  main  sphere. 
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The  "observed"  Ag  values  (due  to  the  disturbance  M)  were  selected  at 
the  surface  of  the  model.  The  first  set  of  observations  (15  observations)  was 
selected  for  to  £ 7'  30",  the  additional  set  of  35  observations  was  selected  for 
8' s oo  £ 25’  (all  35  points  were  located  at  the  surface  of  the  main  sphere).  The 
predictions  of  £ , Ag  and  0 were  carried  out  at  the  surface  of  the  model  (between 
the  observation  points)  for  different  depths  to  the  Bjerhammar  sphere  (Table  5. 1). 
The  RMS  errors  at  the  optimal  depths  are  summarized  in  Table  5.2.  Methods 
1 and  4 give  good  agreements  with  the  theoretical  values.  The  methods  2 and  3 
are  not  suitable  for  this  local  model.  The  importance  of  the  choice  of  radius  of 
the  Bjerhammar  sphere  is  obvious  from  Table  5.1. 


Table  5.  2 

RMS  Prediction  Errors  at  the  Optimal  Depths 
to  the  Bjerhammar  Sphere 


15  observations 

50  observations 

Method1 

Depth*" 

/km/ 

c 

/m/ 

Ag 

/mgal/ 

0" 

Depth" 

/km/ 

c 

/m/ 

Ag 

Angal/ 

0" 

1 

3 

0.04 

22.  7 

0.33 

3 

0.03 

2.27 

0.33 

2 

8-12 

0.37 

27.5 

6.20 

- 

— 

— 

3 

3 

- 

11.9 

34.2 

6.  25 

- 

— 

~ 

3 

4 

3 

0.05 

2.23 

0.32 

3 

0.03 

2.21 

0.33 

1 Norm  used 

2 Difference  of  the  radii  of  main  sphere  and  Bjerhammar  sphere 

3 Not  computed 

No.  of  predicted  points  = 24,  RMS  (anomaly  - mean  value)  = 182  mgal. 
Reference;  Sjoberg  (1975,  Table  19.1) 


5.2.2  Molodenskii’s  Model 


The  surface  of  Molodenskii's  model  is  a cone.  Two  spheres,  whose  centers 
are  located  on  the  axis  of  the  cone,  are  taken  as  the  anomalous  masses  (see  Figure  3). 
The  attractions  of  the  spheres  ma  and  mi  on  the  axis  of  the  cone  4050  m above  the 
reference  plane  are  150  and  100  mgal,  respectively  (ma=«  1.458  x 10XB  kg, 
mi  as  6.334  x 1013  kg).  As  there  is  a disturbing  mass  (mt)  above  the  reference  plane. 


-37- 


it  is  clear  that  this  model  is  much  rougher  than  Bjerhammar's  model.  The 
computations  revealed  that  for  all  methods  (norms)  the  depth  to  the  Bjerhainmar 
sphere  should  be  selected  as  close  to  the  reference  plane  as  possible  (10  meters 
was  used  in  the  computations). 

All  observation  and  prediction  points  were  selected  on  the  surface  of  the 
cone.  A summary  of  the  RMS  prediction  errors  is  given  in  Tables  5.3  and  5.4. 

We  notice  that  the  predictions  are  essentially  better  for  a more  regular  distribution 
of  the  observations  (Table  5.4).  Again,  methods  2 and  3 are  less  favorable  than 
1 and  4.  We  also  notice  that  these  RMS  errors  are  of  a magnitude  larger  than 
those  obtained  for  the  smoother  model  of  Bjerhammar. 

A comparison  between  a minimum  norm  solution  and  the  Dir?  ■ me  hod 
for  Molodenskii's  model  is  given  in  section  8. 


Table  5.  3 


Molodenskii's  Model.  RMS  Prediction  Errors 
at  the  Optimal  Depths  to  the  Bjerhammar  Sphere 


Method1 

Depth 

Am/ 

C 

/m/ 

Ag 

/mgal/ 

0" 

1 

10 

0.67 

25.0 

7.58 

13.78 

2 

200-3500 

0.23 

33.  3 

5.48 

6.25 

3 

10 

13.5 

45.5 

14.24 

28.25 

4 

10 

0.63 

25.0 

7.58 

13.76 

1 Norm  used 

(10  m is  the  minimum  depth  used  in  the  computations.)  No.  of 
observations  = 20  (irregular  distribution),  no.  of  predictions  = 19 
[RMS  (anomaly  - mean)  = 81.7  mgalj.  From  Sjoberg  (1975, 
Table  19.2). 


Table  5.4 

Molodenskii's  Model.  RMS  Prediction  Errors 
at  the  Depth  10  m to  the  Bjerhammar  Sphere 


Method1 

c 

/m/ 

Ag 

/mgal/ 

K" 

8" 

1 

0.09 

15.4 

2.45 

2.70 

3 

5.2 

18.9 

3.11 

5.80 

4 

0.06 

15.4 

2.45 

2.70 

1 Norm  used 

No.  of  observations  = 24  (regular  distribution), 

no.  of  predictions  = 19  [RMS  (anomaly  - mean)  = 81.7  mgal}. 

From  Sjoberg  (1975,  Table  19.3). 
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6.  Integral  Formulas  as  Limiting  Cases 


In  this  section  we  are  going  to  study  the  solutions  by  collocation  and  the 
Dirac  method  for  a continuous  field  of  observations  at  the  surface  of  the  earth. 
Two  problems  are  discussed:  the  uniqueness  of  the  solutions  and  the  existence 
of  the  solutions. 

6.1  The  Uniqueness  of  the  Solutions 

The  non-singular  Dirac  method  was  presented  in  (2.25)  for  a finite  set 
of  observations  (we  exclude  terms  of  degree  less  than  2): 


where 


ri 

Agp  = £ u*k  A (P,k)  ; P = 1,  2,  . . . , 


00 

A (P, k)  = y / (2n+l)c  * (rB/rp)n+a  Pn  (cos  0Pk) 


Using  the  solutions  for  uk*  of  (6. 1)  the  disturbing  potential  (T)  may  be  estimated 
by  formulae  (5.2  ) and  (5.6)  (for  Fik=  f tk) : 


Tp  = rB  > uk*  S (P,k) 


where 


V ^(2n+l)c  n*  / rB  \#+a 

S(P,k)=)  — Pn(cos^Pk) 

^ n-1  \rP/ 


and  the  height  anomaly  is  given  by  { cf.  (5.  7) } : 


Cp  = TP/y 


Let  N go  to  infinity  with  a well-behaving  distribution  of  the  carrier  points.  Then 
we  arrive  at  the  following  integral  formulas  from  (6. 1)  and  (6.2)  in  the  continuous 
case: 


(6.1') 


(P)  = ^ JJ A (p* Q)  u*  (Q)  d ag 


C(P)  = 5T  J7  S(P,Q)u*  (Q)  da, 


In  the  special  case  c*=  2n+  1 (u*  = 4g*)  (6. 1)  and  (6.2')  become  Poisson's  and 
the  generalized  Stokes'  formula,  respectively. 

We  now  study  the  solution  by  collocation  in  the  same  wav.  From  (1.2)  we 
get  the  following  intermediate  step  (for  D=  0): 


and 

(6.2') 


JJ  c (P,Q)  X (Q)dO,  , P=  1,2 N 

and 

Cp=  ^7  JJ  S'(P,Q)  X (Q)  d Oq 

where 

X (Q)  » V Xk  6 (Q-  Qk) 

k=l 

6 = Dirac's  delta  function  {see  (2.24a)} 

Qk  = the  foot  point  at  the  internal  sphere  of  the  normal  ; ..ough  the 
observation  in  P„ 

In  the  limit  N-"00  (with  the  maximum  distance  among  the  observations  approaching  0) 
these  equations  become  the  following  integral  equations: 


(6.3') 

Ag(P)  = 

4n 

JJ  c (P, Q)  X(Q)  dag 

and 

(6.4') 

A r 

c'p»  = ^ 

JJ  S*(P,Q)X(Q)  da, 

Formulae  (6.2  ) and  (6.4  ) are  estimates  of  £ (P).  Furthermore,  the  estimates 
seem  to  differ  for  various  choices  of  c „*  and  cn  of  the  kernel  functions.  However, 
next  we  are  going  to  show  that  these  differences  are  merely  apparent  whenever  the 
solutions  exist. 

First,  we  expand  c (P,  Q)  of  (6.3c)  into  spherical  harmonics  fcf.  formula 

(A.l)  } : 


oo  n 


n = S »=-n 


By  inserting  this  expansion  into  (6.3'  ) we  obtain: 


V B„.  (£)  Y 


(P) 


n=  3 ■ = - n 


Tf 
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where 


B n»  = 


c„ 


4 n(2n+l)  J J \ rq 


jj  v..  (Q)  x (Q,  ** 


In  the  same  way  (6.4  ) and  Ag(P)  may  be  developed  into  the  following  series; 

C(P)>  r.y  Vb.Jj  (£ )'*%..  (P) 


n = 2 B ■=  -n 


and 


($r*- 


<p) 


where 


An,=  ^-JJ  Ag*  (Q)  Y„.(Q)dOq 


From  the  above  expansions  we  obtain  for  en  > 0: 


B A nm 


and 


where 


- IJ  s (r',H>?q)  A%*  (Q)  d Cq 


s<r”^)  = 1 frr  (t;)  p"(cos,M 


n=  a 


S (rp,0pq)  is  the  generalized  Stokes'  function.  Thus  we  have  shown  that 
each  collocation  solution  for  C (for  c„>  0)  in  the  continuous  case  equals  Stokes' 
integral  formula. 
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r ■ — 

The  same  result  may  be  shown  for  Poisson's  and  Vening  Meinesz' 
integrals.  Moreover,  the  same  proof  holds  for  the  non-singular  Dirac  method 
in  the  continuous  case  {formula  (6. 1 ) and  (6.  2 )} . The  uniqueness  of  all  solutions 
may  be  regarded  as  a consequence  of  Stokes'  theorem  (Heiskanen  and  Moritz, 

1967,  p.  17),  which  states  that  a function  harmonic  outside  a closed  surface  S 
is  uniquely  determined  by  its  values  on  S. 

However,  the  conversion  from  one  set  of  c,  (cn*)  to  another  is,  of  course, 
not  valid  unless  the  unknowns  (X  and  u*)  exist  and  the  kernel  functions  are  bov  ad  2d. 
This  question  is  discussed  in  the  next  section. 


6.2  The  Existence  of  the  Solutions 


In  sections  3 - 4 we  have  studied  the  stability  of  the  coefficient  matrix  for 
different  methods,  which  is  of  importance  for  the  practical  application  with  a finite 
set  of  observations.  Now  we  ask  under  which  circumstances  there  exists  a solution 
in  the  continuous  case  (infinite  number  of  observations).  As  will  be  shown,  the 
existence  of  a solution  to  any  of  the  integral  equations  is  considerably  dependent  on 
the  smoothness  of  the  observed  field. 


In  Moritz  (1975),  a proof  of  the  convergence  of  least  squares  interpolation 
is  presented  for  an  element  T of  a certain  Hilbert  space  (with  a given  kernel 
function).  However,  Tscheming  (1977a)  showed  that  the  disturbing  potential  of 
the  earth  (T)  is  not  an  element  of  the  Hilbert  space,  associated  with  the  empirical 
covariance  function.  Subsequently,  Moritz'  proof  is  not  applicable  in  this  case. 


In  our  study  we  start  with  Poisson's  integral  equation  for  the  exterior  of 
a circle  of  radius  rB.  All  observations  are  assumed  to  be  located  on  a circle  of 
radius  r.  Then  we  have  { cf.  formula  (4.4)}  : 


(6.5) 


Ag(0j)  = J k (©j , ete)  Ag*  (6k)  d ek 

-n 


where 

(6.  6) 


k <0J» e*>- 


— V sl"leIn<e’-fl*> 

2 it  L 

n=-oo 


Now  we  assume  that  Ag  may  be  expanded  into  a Fourier  series  and  we  try  a 
corresponding  series  for  Ag*: 


(6.7a)  Ag  (0j)  = £ a.e1*05 

■ ~ — oo 
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and 


(6.7b)  Ag*  (0k)  = £ blell'e« 

1=  -<® 


Furthermore,  we  assume  that  both  the  unknowns  (Ag*)  and  the  observations  (Ag) 
are  uniformly  distributed  with: 


0 


) 


2ttj 

m 


j = 1,  2,  . . . , m 


where  m is  the  number  of  observations  (Dirac  method).  This  implies  that  Ag* 
may  be  written: 


and 
(6.  8) 


Ag* (6)  = ^ Ag*(0k)  6(0-  0k) 

k = l 

a 

Ag(6j)  = £ k (0Jt0k)  Ag*(0k)  , j = 1,  2,  ...,m 

k=  1 


i £0 

Let  us  now  regard  (6.8)  as  a linear  filter  with  each  "tone"  b^  e of  Ag  as  input 
and  a^e1^  y of  Ag  as  output.  Then  we  obtain  from  (6.6)  - (6.8): 


H 


l^a^  J/»_  b£  A o (n|  ln2n  ( j— k)/ a k/» 

' 2tt  L L S 6 G 


k=  1 n = -a= 


£ T jlnlginjaV.  ^ el(t-n)aV. 


- £ I l 


k = 1 


_ W l £ s r*  j /» 


2tt 


m U£,.  (s) 
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• •>* 


where  (see  Appendix,  Corollary  A.  2): 


Hence 


»<..<■>-£  s'""'.  (s'"-®* 

P=  -00 


integer  part  of 


U I 


(6.9)  b£=  2T7ax/(m  u£>.  (s)) 


Moreover  we  notice 

lim  U£,0  (s)  = sl£l  , m - ® 


so  that  for  large  m we  have,  approximately: 


(6.9)  t>£  2 tt  /(m  s l£l) 


We  conclude  that  Ag*  does  not  generally  exist  for  a dense  distribution  of  the  observations, 
and  the  convergence  for  Ag*  as  m — 00  is  very  much  dependent  on  the  behavior  of  a^ 
as  i - ® . 

Let  us  assume  that  the  solution  for  Ag*  exists.  The  predictions  of  Ag  (R,  6) 
are  then  given  by: 


Ag(R,6)=£  T k (0,0k)  Ag*£(6k) 

k=i 


co  a oo 


s r. 


1=-od  k=  \ n=-® 


. m r , lie  , . 

- or  ) b£  e vi,»  (s*.  6) 


2n 


l = -CD 


-46- 


where 


and 


Sr  = rB/R 


v /c  Ox  _ V „ l«+Bh  1«30 
v£,a  (Sr,  W)  - ^ Sr  e 

jC-00 


This  sum  is  given  in  a closed  form  in  the  Appendix,  Proposition  A.  2.  Inserting  (6.9) 
we  finally  obtain: 


(6.10) 


A SB 

Ag  (R,  0)  =y  a£ 

•— j 


(s*,0) 

u£,.(s) 


ae 

e 


Furthermore  it  follows  from  Corollary  A.  3 that; 

it  1 

Vfl,„  "*  Sr  , m - ' 


so  that 


lim  Ag  (R,0) 

a — ♦ co 


£=-°° 


Thus^  Ag  is  convergent  in  thelimit  for  R 2 r . (The  convergence  for  R=  r was  proved 
by  Hormander.  See  Bjerhammar,  1974.)  Subsequently,  we  have  shown  that  although 
the  predictions  are  uniformly  convergent  for  R > r,  the  intermediate  solution  Ag*  does 
not  generally  exist  in  the  continuous  case  (m-®).  This  result  is  of  great  practical 
importance. 

Let  us  now  substitute  (6.  5)  by  the  following  general  integral  equation 
[cf.  formula  (2.  5)] : 


(6.11) 


TT 

Ag  (03)  = J k (0j,ek)  U*(0k)  d 0k 
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where 


2n 


k(6,.e,)=  i.  f /S*  s'-'e"'8--9"! 


To  ensure  that  the  kernel  k is  bounded,  we  suppose  that  the  magnitude  of  c n*  satisfies 


(6.12) 


s1"1  < 1/n  for  n > M 


where  M is  an  integer.  The  solution  of  (6. 11)  becomes: 


u*  (0)  = Y b£ 

£ = -oo 


where  {cf.  (6.9)  and  (6.9  )}  : 


b£  = 2tt  a£/(m/c ult,  (s)) 
or,  approximately  for  large  m: 


(6.13) 


b(»2n  aj/(m/cj*  s,,,,) 


From  (6. 12)  and  (6. 13)  it  is  obvious  that  it  is  generally  not  possible  to  select  a set 
c * such  that  u*  is  convergent  for  arbitrary  a^ . 

Next,  we  proceed  to  the  global  equivalence  of  the  above  derivations.  We 
assume  that  the  gravity  field  of  the  earth  can  be  expanded  into  a series  of  spherical 
harmonics  [Ynit(P)]  for  each  point  (P)  on  the  surface  of  the  earth  (cf.  section  6. 1): 


oo  n 


n=  2 a — ~ n 


Ym  (P) 


where  Y„,  (P)  is  defined  in  the  Appendix,  formula  A.l.  Furthermore,  we  expand 
u+  into  a corresponding  series; 


(6.14) 


CD  T\ 

U*  (Q)  = I I Un"  Y"*  (Q) 

n = a ■=- n 
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where  u nB  are  unknown  coefficients.  By  inserting  these  two  series  into  (6. 1 ),  we 
arrive  at  the  following  identity  for  each  spherical  harmonic; 


■JH. 

2n+l 


un«  = 


A 


nu 


or 

(6.15) 

^ n 

Here  we  have  expanded  Pn  (cos  Hi of  the  kernel  function  A (P,  Q)  according  to  the 
Addition  theorem  (A.l); 


n 

pn  (cos  H)H))=  Y Ynt  (P)  Yn.  (Q) 

b = - n 


and  we  have  also  used  the  orthogonality  property  for  the  spherical  harmonics. 

The  same  technique  may  be  used  to  solve  for  the  intermediate  solution  X (Q) 
in  collocation.  We  insert  the  expansion  for  Ag  (P)  and; 


(6.16) 


X(Q)=£  X„,  Y».  (Q) 

n,  o 


into  (6.3'  ),  and  we  generalize  the  auto-covariance  function  by  substituting  the 
empirical  set  c „ of  (1.5)  by  an  arbitrary  set  of  positive  definite  parameters  c *. 
The  result  is: 


(6.17) 


l l 


C n* 

2n+ 1 


/ rB  \n+a 

Itt)  Y-<p> 


Yn,(Q)Y„v(Q)  d = 


n»  a 


Y„.(P) 


In  this  case  we  can  not  simply  apply  the  orthogonality  property  of  the  spherical 
harmonics  to  identify  the  unknowns,  because  rQ  is  a variable.  However,  if  we 
limit  ourselves  to  the  special  case  rq  = r = constant  (all  observations  on  a 
sphere)  we  obtain: 


/r8  \n+2  c * 

Y 

— A 

/ 2n+l 

^ na 

na 

or 

(6.18) 

x„ 

2n+l 
■ c * 

^ n 

A na 

In  summary, 

the  solutions  for 

u*  and  X are; 

(6.19) 

u*  (Q)  = 

A na 

Y„a  (Q) 

n#  a 

and 

(6.20) 

x <Q)  = y 

n,  a 

ill  A 

c»*  nB 

v r8 

y n+2 

) YnB 

If  we  now  consider  that  the  choice  of  c*  is  restricted  by  the  condition  of  bounded 
kernel  (covariance)  functions,  the  convergence  of  (6.19)  and  (6.20)  is  essentially 
dependent  on  the  smoothness  of  A nB  for  higher  degrees. 

Let  us  for  the  moment  assume  that  (6. 19)  and  (6.20)  converge.  By 
inserting  these  solutions  into  the  integrals  (6.  l")  and  (6.3'  ),  we  obtain  in  both  cases 

<6-21)  Ag(P)  = y £ An.  (^)  Yn.  (P) 

n =2  ■=- n 


This  series  is  identical  with  the  exterior  gravity  anomaly  of  the  earth.  Thus  die 
predictions  converge  to  the  true  values,  whenever  the  unknowns  u*  and  X exist  in 
the  continuous  case  (cf.  section  6.1). 

According  to  the  Hilbert-Schmidt  theorem  (Chambers,  1976,  p.  50),  the 
series  solution  (6.21)  for  Ag  (P)  of  the  left  member  of  the  integral  equation  (6.3'  ) 
is  valid  if  and  only  if: 
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1)  c (P,  Q)  = V cn*  (rB/r)^n*a*  Pn  (cos  d)PQ) 


2>  llrf'  -k  II  Xs(Q)doq<» 

where  r is  the  minimum  radius  to  any  point  at  the  surface  of  the  earth.  Using 
(6.20),  the  condition  2)  may  be  written: 


n=a 


where  we  have  used  the  notation: 


C7n3  is  the  anomaly  degree  variance  at  the  sphere  of  radius  r (cf.  formula  (1.5)). 
1)  and  2 ) we  arrive  at  the  following  inequality  for  the  magnitude  of  c n*  : 


(6.22) 


n^n  o. 


In  the  same  way  we  obtain  the  following  conditions  from  the  integral  equation  (6. 1 
the  solution  (6. 19): 


1)  A (P,  Q)  = £ / (2n+l)c  „*  (ru/rf*3  P„  (cos  ti>pg)  < 00 

n = 0 

and 

2)  m-  . t «*«,  ^ (j-ra)  < - 

n=  0 


and  the  corresponding  magnitude  bounds  for  c*  are  given  by: 
(6.23)  n3ffna<C*(Ii)a(n+a)<  n~3 


From 


) with 
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The  inequalities  (6.22)  and  (6.23)  should  be  regarded  as  necessary  bounds  for  c „* 
in  a convergent  solution  of  collocation  and  the  Dirac  method,  respectively.  Wc  notice 
that  both  (6.22)  and  (6.  23)  are  satisfied  for: 


(6. 23 ) 


< 


n 


-6/2 


However  for  the  gravity  field  of  the  earth,  we  have  according  to  Kaula's  rule  of  thumb 
for  the  variation  of  the  potential  coefficients,  that  an  is  on  the  order  of  n-1  2.  Thus 
we  conclude  that  the  solutions  for  u*  and  X are  not  convergent  in  the  continuous  case. 
The  same  negative  result  is  obtained  if  we  replace  the  observations  Ag  of  (6.1)  and 
(6.3)  by  the  disturbing  potential.  In  this  case,  we  arrive  at  the  same  inequalities 
(6.22)  and  (6.23).  The  degree  variances  of  T (Ctn2  (T))  are  of  nr  .aitude  n 3 according 
to  Kaula's  rule,  and  it  is  clear  that  (6.23  ) is  not  satisfied. 

In  all  the  methods  of  Bjerhammar,  the  idea  is  to  determine  a fictitious  field 
u*  (or  4?*)  °n  the  internal  sphere,  and  then  to  use  u*  (Ag*)  in  the  classical  integral 
equations  for  estimating  geophysical  quantities  such  as  Ag,  T,  £ » and  £.  In  collocation, 
we  are  not  restricted  to  the  determination  of  the  corresponding  intermediate  solution 
X as  was  suggested  previously.  For  example,  the  prediction  of  Agp  may  be  expressed 
directly  as  a linear  combination  of  the  observations  (discrete  case): 

0 

(6.24)  Agp  = Y hp!  Ag! 

i=i 


where  hp!  are  the  weights,  which  are  given  by  the  following  discrete  Wiener-Hopf 
equations: 


(6.25) 


CPK  = hpi  c lk  ; k = 1,  2,  . .. , m 
i=  i 


where  c is  the  covariance  function  of  Ag.  In  the  continuous  case  (6.24)  and  (6.25) 
become: 


(6-24>)  4?  (p)  = ~ Jj  h (P,Q)  Ag  (Q)  d C, 

and 

(6.25')  c(P,Q)=^i-  JJ  h (P,  Q ) c (Q  , Q)  d Oq' 
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If  P is  a point  at  the  surface  of  the  earth,  (6.  25)  has  the  solution: 


h(P,Q)  = 6 (P  - Q) 


where  6 is  Dirac's  delta  function  defined  in  (2.24a).  If  P is  a point  outside  the 
surface  of  the  earth,  we  may  assume  that  the  covariance  function  c is  spatially 
homogeneous  and  isotropic: 


e(P.Q)-X  c n (rB/rP  rq)"+2  pn  (cos  <^q) 
n=2 


In  order  to  solve  for  h we  try  the  following  expansion: 


h (P.Q)  = 7 h „ Pn  (cos  tfpq) 

Li 

n=  2 


where  hn  are  unknown  coefficients  to  be  determined.  Let  us  restrict  ourselves  to 
the  spherical  approximation  of  the  earth,  i.e.  rq  = rq'=  r = constant.  Then  we 
obtain  from  (6.25  ): 


hn  = (2n+l)  (r/rP)n+s 


and  the  solution  for  h becomes  a modified  Poisson's  kernel  function: 


h (P,Q)  = y (2n+l)  (r/rP)n+Spn  (cos  0pq)  = 

n=2 


= L (rp 


JLL 


rp  (rp  + r -2rrpcost/)pq) 


^ - (f:)  - 3 (tt)  — 


and  (6.24  ) becomes  Poisson's  integral  (without  the  terms  of  degrees  less  than  2). 
Thus  we  have  found  that  by  carrying  out  the  collocation  solution  by  solving  the 
Weiner-Hopf  integral  equation,  the  correct  gravity  anomaly  is  recovered  on  and 
outside  the  surface  of  a spherical  earth. 
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In  the  same  way  we  can  solve  for  the  disturbing  potential  in  the  exterior 
of  a sphere.  Let  us  use  the  estimator: 

(6.26a)  T(P)=  JJ  h(P,Q)  Ag(Q)  do, 

where  h (P,Q)  is  given  by: 

(6.26b)  k (P,  Q)  = ~ JJ  h (P,  Q ) c (Q',Q)  d CTq' 

and 

00 

k (P, Q)  = r i> Y (r8a/rp  r)"+3  Ps  ,cos  ^P(5) 

L,  n-1 

n = 2 


c (P, Q)  = Y c„  (ra/r)2(n*a)  Pn  (cos  (dp,) 
n=a 


By  trying  an  expansion  for  h in  (6.26b)  (cf.  the  previous  example),  the  coefficients 
hn  are  easily  identified  and  the  solution  for  h becomes: 


h (P,  Q)  = r S (rp.i/jp,) 


where 

oo  ^ ^ 

S (rp,(M  = £ (-£-)"  Pn(cos*pQ) 

n=a 


Thus  the  solution  for  T /y  from  (6.26a)  with  rp=  r is  Stokes'  formula  (cf.  section  6.1). 


6.3  The  Effect  of  Smoothing  (Noise) 

The  solutions  for  u*  and  X in  the  methods  of  Bjerhammar  and  collocation 
may  be  smooth  ad  by  adding  a noise  covariance  function  d (j,k)  to  the  auto-covariance 
function  c (j,k)  ( fv?e  formula  (1.  2)  and  (2. 16)).  Consider  the  example  of  the  previous 
section  with  a continuous  field  of  observed  Ag  on  a sphere  of  radius  r.  Let  us  replace 
the  covariance  function  of  (6.3'  ) by: 
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(6.27) 


c (P,  Q)  = c (P,  Q)  + d (P,  Q) 


where 


d (P,Q)  = d 6 , d > 0 


6 (l/)pq)  = Dirac's  delta  function,  defined  in  (2.24a).  This  noise  covariance 
function  corresponds  to  pure  white  noise.  By  inserting  (6.  27)  into  (6.  3 '),  we  obtain 
in  accordance  with  (6. 17)  and  (6. 18)  (for  rP  = rq  = r ): 


(6.28) 


X n«  — 


A B»  (r  a/ r)n*3 

d + (rB/r)2(n+a)c  n/(2n+l) 


For  large  degrees  (n)  the  coefficients  Xna  can  be  approximated  by: 


X-.*  Ja„,  (Is) 


r„Y>+a 


Subsequently 


M'l 


x»,  Y„»  (Q) 


is  convergent  whenever  the  spherical  harmonic  expansion  of  Ag  has  the  radius  of 
convergence  r. 

The  prediction  of  new  anomalies  from  the  intermediate  solution  X is  given 
by  (6. 16)  and  (6.  3 ' ).  The  result  is; 


Ag  (P)  = £ X ni  (rB7rrp)n+a  ^ Yn.  (P) 


By  inserting  (6.28)  we  arrive  at  the  following  predictor  and  prediction  error: 


(6.29) 


4 (P)  = y La.r,./r>!('*.‘>  . 

n , (2n+l)d  + c„(rB/r)a(n+a) 


/ rB\"+a 

a-(tJ  y-(P» 
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and 

(6.30) 


Ag(P)  -Ag(P)  = -? 


,7.  d + (rB/r)s(n+a)c»/(2n+l) 


(?) 


Yn.  (P) 


Subsequently,  by  selecting  a sufficiently  small  d > 0 , the  prediction  error  becomes 
arbitrarily  small.  However,  in  practice  the  constant  d has  to  exceed  a certain 
minimum  value  due  to  the  requirement  of  X to  be  within  the  range  of  the  computer. 

In  a similar  way  we  may  approximate  the  solution  for  the  Dirac  method. 

The  smoothing  is  first  of  all  justified  when  the  observations  are  erroneous. 
[Cf.  formula  (1.  2)  and  (1.  3).  ] These  predictions  are  statistical!  tased,  which  we 
conclude  from  the  following  application  of  (1.3).  Let: 


Agp  = c f>  (C  + D)'1  Ag 


be_an  estimator  of  the  anomaly  Agp  , which  is  included  in  the  vector  of  observations 
(Ag).  Furthermore,  Ag  consists  of  a signal  Ag  and  its  error  e : 


Ag  = Ag  + e 


where  Ag  is  the  true  anomaly  and: 


E {*}  = 0 

E { } = 


and  E { e eT } = D 

the  probabilistic  expectation 


Now  it  follows  immediately  that: 


Hence 


E {Agp  } = cp  (C  + Df 1 E {Ag  ) = c p (C  + D)"1  Ag 
E {Agp)  i cp  C'1  Ag  = Agp 


and  we  have  proved  that 


is  biased. 
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In  the  continuous  case  we  may  prove  the  bias  in  the  following  way.  Let 
the  observed  anomaly  field  (Ag)  consist  of  the  signal  Ag  and  the  noise  e : 


Ag  = Ag  + e 


where 


Ag(P)  = £ A„  (^-B)  Y (P; 

nvn 

e(P)  = Y,  €»«  (fT%-(P) 


E {f  } = E {f  } = 0 


The  estimator  of  Ag(P)  when  including  the  noise  in  the  observations  is  in  accordance 
with  (6.  29): 

aa  V cn(rB/r)2(n+2)  /r8\n+s 

Ag(P)  _ ) . + »*  + ^ no  ( P) 

£ (2n+l)  d + c n (rB/r)2^  V 


and  the  bias  becomes 


e fAg(p)} - Ag(P)  = ~y -4^ Y""<p) 

»“  d+(re/r)3(n  a)cn/(2n+l)  V rp' 


The  bias  is  due  to  the  fact  that  Ag  itself  is  not  a random  stochastic  process,  only  its 
errors,  e.  (The  expectation  should  not  be  interchanged  with  the  global  average. ) A 
way  to  diminish  the  bias  is  to  subtract  a low  degree  spherical  harmonic  reference 
expansion  from  the  observations  prior  to  the  prediction.  Even  though  we  cannot 
correct  for  all  the  bias,  the  resulting  residuals  are,  hopefully,  more  random  than 
the  original  observations.  The  general  prediction  formula  (1.2)  should  therefore  be 
modified  in  the  following  way  [cf.  Sjoberg,  1975,  formula  (16.1)  and  Koch,  1977, 
formula  (18)]: 


Vi  = V,  + C , (C  + D)-1  (Ag  - Ag) 


(6.31) 


where 


v . , Ag  = low  degree  spherical  harmonic  expansions  for  vt  and  A g. 


This  is  the  well-known  collocation  with  the  inclusion  of  systematic  terms  represented 
by  v , and  Ag . We  conclude  that  this  formula  should  be  used  not  only  in  local  studies, 
but  also  in  global  applications  of  collocation  to  reduce  the  bias  in  the  model. 


6.  4 The  Limiting  Case  of  a Least  Squares  Solution 

Finally,  we  are  going  to  study  the  convergence  of  the  least  squares  solution 
of  Poisson's  equation  for  the  circle.  We  assume  that  we  are  using  the  Hirac  method 
with  m>N.  The  residuals  may  be  written: 

N 

€j  = Ag j - £ k (6j,ek)  Agk*  , j = 1,  2,  ....  m 

k = 1 


The  minimum  of  the  square  sum  of  the  residuals  is  obtained  for  the  normal  equations: 

D 0 N 

y k (0  j ,e,)  Ag  J = y £ k (0  j , 0q)  k (0, , 0k)  Ag  f 

J=1  J=l k=l 

where 

q = 1,  2,  ...,  N 

Insc  rting 

00 

O 0 

an  e1"  J , 6j  = 2TTj/m=  2TTr  j/N 

n=  - oo 

and 

00 

Ag*  (6*)  = y bj  eaQ*  , ek=  2 nk/N 

£=-°° 

where  r is  an  integer  (m  is  a whole  multiple  of  N)  into  the  normal  equations  it  can 
be  shown  that: 

b£=a*y  s,e+Npi/Ny  y ^ 

P P q 
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r 


where  all  summation  intervals  are  [-«,  +«].  Specially,  for  m - » we  arrive  at; 

y S'*+Np' 

(6.27)  a* 


"I 


3 IP.+  Npl 


where  the  sums  are  given  in  a closed  form  in  the  Appendix,  Corollary  A.  2.  it 
follows  from  (6.32)  that  Ag*  will  generally  not  exist  when  N approaches  infinity, 
because  for  large  N we  have  approximately: 


s*  a£  s l{,/N 


The  prediction  of  an  arbitrary  point  Ag  ( 0 , R)  from  Ag*  and  (6.32)  gives; 


r ..it  ♦ %£»! 


Ag  (fl.R)  = Y a,  6 V,.,  (6.  R)  £ 


T-*"* 


Npl 


where 


^ IMDI  1 Nn  V 

Sr  e , Sr  = re/R 


Specially  for  8=  9 } = 2n  j/N  and  R = r we  obtain: 


(y  a11***) 

A .iQ  '4*  ' 

Ag  ( 6j , r)  = ) aAelWJ 


y S2,(Un' 


and  finally 


y-C  0 

_ a£  e1  3 = Ag  (8}  , r) 


Hence,  the  least  squares  prediction  converges  to  the  true  value  when  N approaches 
infinity. 
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7.  On  the  Choice  of  Radius 


In  the  studies  of  the  Bjerhammar  problem  little  attention  has  been  paid 
to  the  choice  of  radius  of  the  Bjerhammar  sphere  (rB).  However,  in  Sjoberg  (197b) 
the  importance  of  this  choice  is  demonstrated  for  different  covariance  functions. 

In  this  section  we  are  going  to  use  the  following  one -dimensional  prediction  estimate 
on  a circle  as  a basis  (cf.  the  previous  section): 


(7.1) 


where 


4?  (r,0) 


vi» » 

u£.« 


i£  0 

e 


V£,  a — V£,  u (S 


,e>  = f 


J=-00 


U£,b  = V£,  B (S,0) 


and  m is  the  number  of  observations.  The  corresponding  prediction  errors  are; 


(7.2a) 


wher; 


e(0)  = Ag(r,  0)  - Ag(r,  0)  = Y e,(0) 


£=-« 


(7.2b) 


1 16 


We  define  the  mean  prediction  error  in  the  following  way: 

2tt 

(7.3)  ? = ~ J €<®>  d 6 

0 


First  we  prove  the  following  proposition. 


-60- 
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Proposition  7.1:  If  a 0 = 0 then 


oo 

< |Y  a 


0 < |e  | < I ) a, 
1=  —00 


and 


f - 0 

e-V  ajB  , 


s - 1 
s - 0 


where  m is  the  number  of  observations. 
Proof:  Using 

2n 


~ J Ag  (6)  d 6 = a0  = 0 

o 


it  follows  from  (7.1): 


2n  2n 

= — r 4(e)d9=y  is. y sb*-d  i [,<■•■-)«« 

2 TT  J 7*  U^*  “ T1  2 J 


V — = Y ajB 

— > U Jo  » a 1+  ® — 1 

J J 


Here  we  have  used  Corollary  A. 2 with  1 = jm 


■*Js  • a 


1+S» 

1-s" 


From  this  expression,  the  proposition  readily  follows. 

Although  s-1  gives  the  mean  prediction  error  0 and  has  also  proved  to 
give  the  most  stable  solutions  (section  3),  the  "best"  choice  of  radius  tb  must  be 
determined  in  some  mean  square  sense.  Let  us  therefore  study  the  prediction  of 
the  signal: 
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* 


(7.4) 


Ag£(0)  = b cos  (£0)  + c sin  (£0) 


r 


at  the  circle  of  radius  r.  We  define  the  relative  variance  in  the  following  wa> 

2tt  2tt 


(7.5) 


where 


Re  = J f2  (0)  d 0/  J Ag/  d 0 


e(0)=  Ag£<0)  - Ag^(9)  (prediction  error) 


It  is  shown  in  the  Appendix,  Proposition  A.  4,  that  (7.5)  may  h ■ written  'for  £>0): 


(7.6) 


where 


(7.6  ) 


R* 


if  £ / yP,  p = 1,  2,  3, 


1 + 2 0 1 

— --—  + 


S;  if  £=  — P 


2(1+0)  2(1  + 0) 


So  = 1 + 


l-s”  1+s 


2 m - 4 n 


1+s"  (l+sB_an)a  l + S*-a" 


1-S-  2,i-. 


[-1 


n = i ' | — 1 m 
. m 


0 = (c/b)‘ 


From  (7.6)  and  (7.6  ) we  draw  the  following  conclusions; 


a)  lim  Re  = 0 , m -* a 

b)  lim  S£  = 1 , s - 1 


c) 

d) 

e) 
0 

g) 

h) 


f £ < m/2  then  lim  R£=0,  s-0 

f £ = m/2  then  S £>  1/2  and  lim  S^=  1/2,  s - 0 

f £ > m/2  then  Se>  1/2 

f £ = mp  then  S£<2  and  lim  S^=  2,  s-0 

f £ = m(p+  1/2)  then  S£<  3/2  and  limS£=3/2,  s-0 

f £ > m/2  and  £ / mp/2  then  lim  = 00 , s - 0 
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RELATIVE  VARIANCE  (%) 


/ /=50(=m) 


* = 30 


- / / 


i=25(  = m/2) 


V.V.  i = l5 

x.  • • . I ~ • • • 


The  relative  variance  Sx  in  the  Dirac  method  for  the  circle,  using 
Poisson's  kernel,  is  given  for  signals  of  various  frequencies  (1)  as 
a function  of  h/a,  where  h=  distance  from  the  Bjerhammar  circle  to 
the  circle  of  observations  and  a*  spacing  of  the  data  (2  nr/m).  m=50. 


(h/a)  optimum 


h- 

0.05  0 


Relative  variance 


The  optimum  ratio  h/a  and  the  corresponding  relative  variance 
(Ri)  given  as  functions  of  l/m,  when*  1 is  the  frequence  of  the 
signal  and  m is  the  number  of  observations  (1?  0.5  m). 


These  statements  are  illustrated  in  Figures  4-5.  It  is  obvious  that  fur 


low-frequency  observations  (X  £ m/2),  the  radius  rB  should  be  chosen  as  small  as 
possible  in  accordance  with  the  numerical  difficulties  discussed  in  sections  3-4  . 

(See  also  section  8.)  For  higher  frequencies  (X  > m/2),  the  relative  variance  it£ 
is  at  least  50%  (see  also  Bjerhammar,  1977a).  Accordingly,  signals  of  frequency 
higher  than  m/2  can  not  be  predicted  in  a satisfactory  way. 

In  a real  prediction  case  we  do  not  have  the  above  ideal  situation  w :t,h  a 
uniform  distribution  of  data  of  one  single  frequency,  and  the  most  favorable  distance 
rB  is  a function  of  the  unknown  coefficients  (a^)  and  of  the  mean  spacing  of  the  dat  . 
In  most  cases  empirically  derived  covariance  functions  are  associated  with  a specifi 
radius  of  the  Bjerhammar  sphere.  Several  tests  in  Sjoberg  (1975)  revealed  that  the 
solutions  by  collocation  are  rather  unsensitive  to  the  choice  of  kernel  functir  m , if 
the  optimal  depth  to  the  sphere  is  chosen  for  each  of  them.  Hov  er,  the  choice  of 
radius  (changing  for  each  distribution  of  the  observations)  is  mportant  for  a good 
result.  In  conclusion,  a study  of  the  relation  between  the  optimum  radius  rs  and  the 
spacing  of  the  observations  is  recommended  for  each  specific  kernel  function.  For 
instance,  in  the  numerical  example  of  the  next  section,  the  optimum  depth  to  the 
sphere  in  the  Dirac  method  is  approximately  half  the  distance  between  neighboring 
observations. 


8.  Computations 

The  iterative  method  described  in  section  4 will  be  used  to  demonstrate  the 
rate  of  convergence  of  the  solutions  for  Ag*  and  X in  the  Dirac  method  and  collocation, 
respectively  [see  formulae  (3. 1),  (3.  2)  and  (4.  2)].  The  stability  of  each  method  is 
reflected  by  the  number  of  necessary  iterations  for  the  solution  of  the  vector  Ag*  or  X. 

First,  the  case  with  the  degree  variances  equal  to  2n+  1 will  be  considered. 
The  observations  consist  of  87  free  air  gravity  anomalies  regularly  distributed  with  a 
spacing  of  approximately  0°.  5 (Table  8.1).  The  50  prediction  points  and  anomalies 
are  given  in  Table  8.2.  The  moan  value  (Ag)  of  these  50  anomalies  is  -2.8  mgal  and 
the  [IMS  value  of  Ag  - Ag  is  ±13.4  mgal.  In  the  computations  we  assume  that  the 
mean  sea  level  radius  is  6370  km.  The  prediction  results  are  given  in  Table  8.3  and 
Figure  6.  From  the  table  we  notice  that  the  number  of  necessary  iterations  increases 
with  the  depth  to  the  Bjerhammar  sphere  and  that  this  increase  is  more  pronounced  for 
collocation  than  for  the  Dirac  method.  The  solution  by  the  Dirac  method  for  one  specific- 
depth  is  identical  with  the  collocation  solution  for  half  that  depth.  This  result  implies 
that  the  latter  method  is  twice  as  sensitive  to  changes  of  the  depth  as  the  former. 

For  comparison  we  study  also  the  predictions  using  collocation  with  the 
empirical  covariance  function  implied  by  subroutine  COVA  (Tscherning  and  Rapp,  1974). 
The  computations  are  performed  for  different  lower  bounds  (N,ln)  of  the  covariance 
function.  In  each  case  the  RMS  prediction  error  10. 1 mgal  was  obtained.  The  number 
of  necessary  iterations  are  given  in  Table  8.4. 
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Table  8.2 


Prediction  points.  50  free  air  gravity  anomaly 
stations  in  Manitoba,  Canada. 


HO 

LOHCITUDE 

LATITUDE 

ALTITUDE 

AHOP1ALY 

DEG 

DEG 

METER 

MGAL 

16009 

50.59499 

91 . 14033 

386. 181 

° . 76 

1601  1 

50.74666 

90.96165 

379 . 476 

>.63 

10499 

50.79678 

92.36719 

425.501 

8.82 

16042 

50.95000 

9 1 . 75833 

402.031 

-9.20 

16012 

50.88333 

91 .01666 

378.257 

1.71 

10052 

50.59193 

92.61716 

428.854 

12.66 

10047 

50.22284 

92.84720 

358.445 

-3.57 

10030 

50.21497 

92.37631 

357.833 

-3.37 

15372 

50.49666 

91.87166 

357.225 

12.68 

273 

50.40833 

9 1 . 50833 

370.027 

3.99 

15370 

50.75499 

9 1 . 88333 

391.973 

7.53 

'6047 

50 . 68333 

91.41998 

382.219 

13.54 

10496 

50.80524 

92. 14000 

396 . 545 

-2.  18 

100BI 

50.21280 

93.23965 

364.845 

7.53 

625 

50.30666 

93. 17999 

361 .493 

4.41 

16002 

51 . 12833 

90 . 86833 

380.390 

-21.34 

16050 

50.29498 

91 . 38998 

373.685 

- 14.68 

15331 

50.74500 

90.75665 

390.449 

7.83 

15334 

50.42332 

90.71666 

391.058 

2.09 

16056 

50. 18333 

90.68832 

404 . 774 

-15.16 

15384 

50. 12999 

90.23999 

417.576 

-3.35 

9 164 

48.91499 

90.60001 

457 . 809 

-11.23 

13576 

49.02196 

91.96181 

423.501 

-21.95 

15715 

49 . 84846 

90.40294 

442.265 

-2.80 

15723 

49.74696 

90.75410 

435.559 

-5.  18 

15675 

49.86766 

92.09277 

382 . 524 

-23.51 

10019 

49 . 84630 

92.39220 

366.674 

-18.90 

15378 

49.92999 

9 1 . 38333 

388.925 

-29.92 

10203 

49.76701 

94.87744 

359.969 

2.  1 1 

1 1 455 

49.62477 

94.02734 

373.990 

-10.26 

10222 

49.59979 

94.35619 

323 . 088 

5.44 

5710 

49.43166 

96 . 27499 

357.225 

5.27 

5546 

49.72333 

95.24666 

338.328 

16.62 

10198 

49 . 62842 

95 . 49942 

318.211 

-1.36 

5531 

49.71t>66 

94.93666 

359 . 664 

10.92 

5076 

49.71500 

94.80666 

345 . 643 

14.  12 

10078 

49 . 90807 

93. 14622 

372.770 

-8.51 

10104 

49.68629 

93.87337 

379. 171 

-18.87 

10098 

49 . 30528 

93.51170 

338.023 

6.60 

5084 

49.81 332 

92.97501 

355 . 092 

-10.  10 

10014 

49.62167 

92.44067 

385 . 572 

18.04 

12002 

49.49670 

92.69481 

403.555 

11.76 

10007 

49.22151 

92 . 46306 

405 . 384 

-23.49 

213 

49. 14166 

92.70332 

388.620 

-22.67 

15725 

49 . 70807 

9 1.09705 

435 . 559 

-8.48 

15736 

49.66466 

91.81667 

419.405 

-24.77 

15738 

49 . 24759 

9 1 . 46384 

445.617 

-13.99 

10249 

40.66917 

93 . 27843 

337.718 

-6.09 

13652 

48.56082 

90.73390 

447.751 

-23.06 

15748 

48.83151 

90.96706 

445.922 

-18.72 
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RMS  prediction  errors  compared  for  various  depths  to  the  Bjerhammar  sphere. 


DEPTH  TO  THE  BJERHAMMAR  SPHERE  (km) 


Table  8.  3 


Comparison  Between  the  Predictions  of  50  Free  Air 
Gravity  Anomalies  From  87  Observed  Anomalies 
Using  the  Dirac  Method  and  Collocation* 


Dirac 

Method 

Collocation 

ra 

(km) 

depth 

(km) 

No.  of 
iterations 

RMS  error 
[mgal] 

No.  of 
iterations 

RMS  error 
[mgal] 

6315 

55 

30 

10.7 

30 

11.8 

6320 

50 

25 

10.6 

?n 

11.7 

6325 

45 

16 

10.5 

30 

11.5 

6330 

40 

10 

10.4 

30 

11.3 

6335 

35 

7 

10.2 

30 

11.1 

6340 

30 

5 

10.2 

30 

10.9 

6345 

25 

4 

10.2 

25 

10.6 

6350 

20 

4 

10.4 

11 

10.4 

6355 

15 

3 

11.0 

5 

10.2 

6360 

10 

2 

12.0 

4 

10.4 

6365 

5 

2 

13.2 

2 

11.9 

6370 

0 

l 

13.5 

i 

13.5 

* The  iterations  of  Ag*  and  X are  interrupted  whenever  the  RMS  residuals 
are  less  than  0.25  mgal. , the  maximum-residuals  are  less  than  0.5  mgal. 
or  the  number  of  iterations  exceed  30.  Degree  variances  = 2n+l. 


Table  8.4 


Number  of  Necessary  Iterations  of  the  Vector  C-1  Ag 
to  Satisfy  Either  of  the  Conditions  RMS  Residual  < 0.25  mgal 
or  | max.  Residual  |<  0.5  mgal.* 


N.in  3 5 3 13  20 

No.  of  iterations  18  16  14  13  10 


* Subroutine  COVA  is  used  for  different  minimum  degrees 

(N  ain  ) 
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The  empicical  covariance  function  in  collocation  gives  a slightly  better 
prediction  result  than  the  Dirac  method  for  its  optimum  radius  of  rB.  This  gain 
is  obtained  at  the  cost  of  at  least  twice  as  many  iterations.  For  a denser  spacing 
of  the  data,  the  number  of  necessary  iterations  will  increase  and  the  difference 
between  the  two  methods  becomes  more  pronounced.  With  this  reasoning,  the 
non-iterative  solution  may  fail,  when  using  collocation  (due  to  numerical 
singularity),  while  a solution  by  the  Dirac  method  may  still  be  useful.  For  example, 
the  prediction  result  reported  in  Table  8.5  shows  small  differences  between  collocation 
and  the  Dirac  method  (93  well  distributed  observations).  However,  when  a few  more 
observation  points  were  included  (99  points)  between  the  previous  ones  (with  minimum 
spacing  675  meters),  the  RMS  prediction  error  of  £ increased  to  11 '.'4  for  the 
collocation  type  of  solution,  while  the  RMS  error  of  the  Dirac  method  was  still 
useful  (1!'3).  (From  Sjoberg,  1975.) 


Table  6,  5 

RMS  Prediction  Errors  for  Molodenskii's  Model 


No.  of 

Collocation 

Dirac  Method 

obs. 

C 

Ag 

C Ag 

i" 

/m/ 

/mgal/ 

/m/  /mgal/ 

93 

99 

0.01 

1.4 

EBB 

0.12  1.8 

0.90 

1.3 

* Depth  to  the  Bjerhammar  Sphere  = 10m.  No.  of  prediction  points  = 37. 
c„*=  2n+l.  Reference:  Sjoberg  (1975,  section  19.  2.  3). 

9.  Conclusions 

The  purpose  of  this  report  has  been  to  compare  some  methods  of 
A.  Bjerhammar  with  collocation  for  the  solution  of  the  boundary  value  problem  in 
physical  geodesy.  In  practice,  the  number  of  observations  are  finite  ("discrete 
boundary  value  problem").  In  the  present  application  collocation  is  most  frequently 
identified  as  Wiener-Hopf  prediction  of  stochastic  processes,  in  which  case  the 
covariance  functions  are  assumed  to  be  known  (homogeneous  and  isotropic).  The 
main  problem  is  therefore  to  find  the  appropriate  covariance  functions.  Rigorously 
this  has  been  proved  to  be  an  impossible  task,  because  of  the  non-ergodicity  of  the 
empirical  covariance  functions  (Lauritzen,  1973). 
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Stokes'  formula.  For  a finite  set  of  observations  there  is  always  a fictitious  field 
(Ag*)  at  an  internal  sphere  that  satisfies  the  integral  equation.  In  Bjerhammar's 
applications  Ag  has  originally  been  considered  as  mean  anomalies  over  blocks  of 
certain  sizes  at  the  internal  sphere  (the  Bjerhammar  sphere).  Dependent  on  the 
number  of  such  blocks,  different  types  of  solutions  of  the  problem  are  obtained. 

If  the  number  of  blocks  (N)  are  less  than  the  number  of  observations  (m),  a unique, 
solution  is  given  from  adjustment  by  elements.  For  N > m a unique  solution  is 
given  by  condition  adjustment,  which  solution  minimizes  the  norm  of  the  unknowns 
(Ag*).  In  the  special  case  N - 00  and  well-behaving  surface  elements  we  have  shown 
that  the  minimum  norm  solution  of  Bjerhammar  ( , Ag*  ||=  minimum)  approaches  the 
solution  by  collocation  with  the  degree  variances  equal  to  2n+l.  Furthermore,  this 
proof  has  been  generalized  to  yield,  that  for  each  set  c n in  collocation,  there  is  a 
corresponding  minimum  norm  ||u*  n in  the  generalized  Bjerhamn  .r  approach  (section 
2.2).  Thus  the  problem  of  selecting  the  degree  variances  in  collocation  is  identical 
with  the  problem  of  selecting  the  minimum  norm  in  the  Bjerhammar  method.  It 
should  be  stated  that  already  Krarup  (1969)  regarded  collocation  as  a generalized 
approximation  for  a specified  minimum  norm. 

A different  type  of  solution,  called  reflexive  prediction,  was  introduced 
by  Bjerhammar  (1974).  in  this  method  the  external  gravity  field  is  assumed  to  be 
generated  at  a priori  selected  fictitious  carrier  points,  on  or  outside  the  Bjerhammar 
sphere.  Once  the  carrier  points  and  the  sphere  are  defined,  Poisson's  integral 
equation  can  be  applied  rigorously,  and  the  result  is  a set  of  linear  equations.  If  all 
carrier  points  are  located  at  the  internal  sphere,  the  method  is  called  the  Dirac 
approach.  Due  to  the  arbitrary  location  and  number  of  carrier  points,  a wide  variety 
of  solutions  are  possible.  Of  special  interest  are  filtering  (less  number  of  carrier 
points  than  observations)  and  the  non-singular  Dirac  method  (see  below). 

In  this  report  we  have,  first  of  all,  compared  collocation  and  the  non- 
singular Dirac  approach  with  carrier  points  located  at  the  intersections  of  the  internal 
sphere  with  the  radius  vectors  of  the  observations.  The  coefficient  matrix  of  the 
latter  method  is  non-symmetric  in  contrast  to  the  former.  Moreover  it  was  shown 
in  section  2 that  for  a constant  geocentric  radius  of  all  observations  in  an  area  (and 
cn  = 2n+  1)  the  two  methods  give  identical  prediction  results  for: 


(9.1) 


h0  = 2 he 


where  h0  and  he  are  the  depths  to  the  Bjerhammar  sphere  in  the  Dirac  method  and 
collocation,  respectively.  This  result  has  been  verified  to  hold  also  for  approximately 
constant  geocentric  radius  of  the  observations  (section  8). 
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It  was  demonstrated  in  sections  3 and  4 that  the  numerical  stability  of 
the  solutions  differs  for  the  two  methods  (collocation  and  Dirac).  For  a given 
radius  of  the  Bjerhammar  sphere  the  Dirac  method  is  more  stable.  In  the  numerical 
example  of  section  8 it  was  found  that  for  Poisson's  kernel  (c  * = 2n+l),  the  con- 
ditionings of  the  two  methods  are  equal  provided  that  formula  (9. 1)  is  satisfied.  For 
this  kernel  function  the  solution  by  collocation  is  twice  as  sensitive  to  the  choice  of 
Bjerhammar  sphere  as  the  solution  by  the  impulse  method,  a fact  that  is  important 
for  the  solution  of  a large  system  with  a dense  distribution  of  the  observations. 

In  the  continuous  case  (with  observations  covering  all  the  surface  of  the 
earth)  we  have  found  (section  6)  that  all  predictions  with  the  generalized  Dirac 
method  and  collocation  are  unique  for  various  sets  of  positive  c if  (cn),  whenever 
the  solutions  exist.  This  result  may  be  regarded  as  a consequence  of  Stokes' 
theorem.  In  general,  however,  the  intermediate  solutions  u*  and  X do  not  exist 
in  the  continuous  case.  The  existence  of  the  solutions  requii*es  that  the  degree 
variances  of  the  observations  are  at  most  of  magnitude  n~°,  a condition  which  is  not 
satisfied  for  the  gravity  field  of  the  earth  according  to  Kaula's  rule.  Approximate 
solutions  may  be  found,  for  example,  by  adding  a positive  constant  to  the  kernel 
• function.  In  the  same  way  it  was  found  in  section  6.3  that  a solution  may  exist  when 

considering  that  the  observations  are  erroneous.  However,  these  solutions  are 
statistically  biased.  If  collocation  is  carried  out  by  solving  the  Wiener-Hopf  integral 
equation,  a convergent  solution  is  obtained  outside  a sphere  (all  observations  on  the 
sphere).  However,  inside  the  bounding  sphere  of  the  real  earth,  the  convergence  is 
still  not  proved. 

It  should  be  noted  that  the  deterministic  approaches  by  Bjerhammar  through 
Poisson's  and  Stokes'  formulae  do  not  provide  estimated  prediction  errors,  as  is  the 
case  in  the  stochastic  process  approach  (i.  e.  collocation  according  to  Moritz).  How- 
ever, as  the  uncertainty  of  the  covariance  functions  has  a direct  impact  on  the  error 
estimates,  these  estimates  might  be  of  limited  value. 

For  large  systems  the  free  choice  of  carrier  points  in  reflexive  prediction 
might  be  advantageous  (filtering  is  possible).  However,  experiences  by  Bjerhammar 
(1977b)  indicate  that  numerical  difficulties  may  occur  in  the  filtering  process,  due  to 
ill-conditioning. 

Finally,  we  like  to  mention  that  the  original  approaches  of  Bjerhammar 
are  designed  to  solve  problems  in  physical  geodesy,  first  of  all  the  geodetic  boundary 
problem,  for  gravity  anomalies  as  the  only  source  of  data.  However,  it  is  not  difficult 
to  modify  the  methods  in  order  to  take  other  types  of  geophysical  information  into 
account.  Thus,  both  collocation  and  Bjerhammar's  methods  possess  a flexibility  for 
the  processing  of  heterogeneous  data. 
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10.  Extensions  and  Recommendations 


This  study  has  shown  that  reflexive  prediction  (the  Dirac  method)  is  less 
sensitive  to  changes  of  the  radius  of  the  internal  sphere  than  collocation  for  the 
particular  kernel  and  covariance  functions  with  c * = c n = 2n+l.  The  prediction 
results  and  the  stability  of  the  two  methods  are  the  same,  if  collocation  is  applied 
with  half  the  depth  to  the  Bjerhammar  sphere  used  in  reflexive  prediction.  We 
recommend  these  comparisons  to  be  carried  out  for  other  covariance  functions. 

Of  special  interest  would  be  to  compare  an  empirical  covariance  function  with  a 
best  fitting  kernel  function  (cf.  Sjoberg,  1975,  section  18). 

Special  attention  should  be  paid  to  the  relation  between  the  distribution  of 
the  observations  and  the  optimum  radius  of  the  Bjerhammar  sphere  for  /arious 
covariance  (kernel)  functions. 

A procedure  to  estimate  the  prediction  errors  in  reflexive  prediction  is 
of  interest  for  the  user  of  the  method.  In  order  to  reduce  the  bias  of  the  predictions, 
we  recommend  the  use  of  formula  (6.31)  in  all  applications  of  collocation  to  the 
geodetic  boundary  problem  and  related  problems  in  physical  geodesy.  A corresponding 
formula  should  be  used  in  reflexive  prediction,  if  the  kernel  function  is  modified  to 
take  noise  into  account. 

Theoretically,  it  is  of  interest  to  reveal  whether  the  Wiener-Hopf  type  of 
predictions  [formulae  (6.24'),  (6.25'),  (6.26  a-b)]  converges  when  applied  to  a 
continuous  field  of  observations  at  the  surface  of  the  earth. 
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(A . 2)  ~ ^ r, ■ Y n i d O 6 ,n  8 i| 

a = the  unit  sphere 
5 = Kronecker's  delta 


Corollary  A.  1; 

(A.  3)  -—  | Pn  (costt  iic)  P3'  (cos  J)()  d crk  = — P„(cos  $ „)  x 6nn' 

4rt  J J Jn+1 

The  corollary  follows  from  (A.  1'  and  (A.  2). 

Proposition  A.  1 : 

If  A is  a matrix  of  dimensions  (m  x N),  where 

00  n+3 

(A)Jk  = -^-r  /(2n+l)  -/cT  (— ')  Pn  (cos0Jk)  Aak 
4n  j v r3  / 

n=o 
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and 


and 

N 

V AOn  = 4 n 

K=  1 

then 

00  2 n^*  3 

lim(AQAT)u  = Y cn  ( ^ Pn(cos^u) 
->  ' r,rj/ 

N — n = 0 

na  x ACT”4  O 
k k 


P roof : 


N 

(AQAT)U  = 4n  ^ (A)ik  (A)Jk  Aok_1  = 

k=i 


N co  n + 3 00  n+ 2 

-±y  y /(2n+l)  Ac „ (12)  P.(cob^u)V^(2S5)^7(^)  pn' (cos  0Jlt)  Aot 


k = l n=0 


n = 0 


For  N in  such  a way  that  max  Aak  - 0 the  summation  over  k becomes  the 
corresponding  integral,  so  that: 


lim  (AQA1),.,  = -2-  V V v/(2n+l)(2n+l)  (•?  ) (t5)  x 

4n  j j v r */  \ r j / 

n - ® n n' 


x J J pn  (COS^k)  P„-  (COS^Ijk)  d CT* 
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Corollary  A.  1 finally  yields; 


lim  (A  Q At)u  - Y c„  Pn(cos^u) 

M —oo  -*  ' TlTj  / 


Proposition  A.  2 


(A. 4)  V (aeI(fl  ) M-tt  N 


, N-ai«i*-  zn'y  _in <n  _ 

1 i £ e ~ 

1_sNew<a  l-sNe'1N<£)  J 


for  0 s s < 1 and  a = I = integer  part  of  ^ and  <a  - { ® if  m > 0 

L N J N V 1-6  ifiTK'  0 


Proof: 


Let  us  substitute  j of  (A. 4)  by  k - Ct.  Then  we  obtain  in  the  left  member; 


^ gla  + Nk-NCUl  el(«  + M(-N0f)t)  g + g 


where 


00 

Si  = (se^M-a«  V (se‘Or  = (se'0)|'|-<vN  / [ i-(Se^)N  ] 


CO 

Sa  = (se“  1 )N  Of  — |o,  | y(se-^r  = (se-^a-  u I /[i_(se-«^N 


Corollary  A.  2; 


(A.  5) 


s ’ +NJ  L (S  '*  '”®N 


) / (l~sN  ) 


The  proof  follows  directly  from  the  proposition  for  f)  = 0. 
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Corollary  A.  3, 

(A.  6) 


lim  7 s'»+nji  e‘(»*^)<p  = s,.,ei^5 

N -*oo  j—  -00 


lhe  proof  follows  directly  from  (A. 4)  when  considering  that  a-O,  N- 


Proposition  A.  3; 

If 


AJk  = 1 Y s|n|eln<6J-  6k> 
J N L 


and 


then 


= 2tt  k/N 
0 j = 2 n j /m 
N = pm  p = integer 

= 0 


m n \ 

S (m,  N,  s)  = N V Y Aq„  AJk  - 2N  V A 


J=1  k=l 


k = 1 


Qk 


equals 


S(m,N,s)  = r(2-m)(l+sN)  + 2m  — - 

*—  / -*  _■ 


-s2N-S"  i_s2N-s 

- 4 


(1-S2*)(1-S3N")  (l-s2)(l-s2N'1) 


+ N 

1-s*  1-S  J 


Proof: 


For  0,  = 0 we  have: 


-78- 


(A. 7)  Bj  = NSY  Aqk  Ajk  = V V s'n'  + ,eielfejt  e~‘(n  + 

k=  1 n=-oo  £=-oo  )<—  i 


= nY  Y s|I,+  *n+Nrl  e1(n  + Nr)0.  = NV  s.n.  el-6-  y Sm  + Nr| 


n=-cor=-oo 


n = -oo  r=  — oo 


and 


(A. 8) 


CD  CD 


Y Bj  = N mV  s-Hy  sMt  + l’ 

L-*  —i  -i 

J=1  t =-00  r=-C0 


Inserting  (A.  5)  into  (A.  7)  for  0j  = 0,  we  obtain: 


00  00 


nY  Alk2=Y  s|nlY  s|,1+Nr,=  -1-  Y (s2,n,“aN  + sN  + ®\  = 


k = 1 n = — oo  r=  — oo 


1-S  -L 


00  Nk+N-1 


[-i + 2y  y (s2—+  = 


k=0  n=Nk 


1 r-  -N.oV 


-1 -s  + 2 


rx_o  “ i 

k=o 


+ 2 N Y s' 

k — 0 


N + 


= ! r.l  V+  Zt1-8-  > + 2Ns"  1 

1-s”  L + ~ J 

In  the  same  way  we  obtain  from  (A. 8); 

— V n - m r,  „N  ,o  l-sa"  -s‘  . 2psN  o 

N l B‘  i-s"  L 1 ■,+2,T?W-)  "T?] 
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so  that 


S (m,N,s)  = -J-  V Bj-  2nY  A„f 

N j -j 


= f-m  -msN 


+2m-Jj:^ — i^r+  + 2 + 2sn  - — 

(l-S  )(1-S3N  ) 1-8*  (l-S^l- 


-J(i-g3W~3)  _ 1NsnT  /n_s^ 

(l-S^l-S3"'1)  1-8  j 1 


Proposition  A. 4 


where 


Ag  ( 0 ) = b cos  (l  0 ) + c s in  (l  0 ) , l > 0 


^(6)  = IT-  e'£6+  — e_li,e 

u£»  ■ u£, . 


a£  = i(b  - ic) 


a-  ( = ^(b+  ic) 


where 


(0)»  T s'"j+*vj6 


U£,  . = Vh  . (0) 

3TT 

l(A{i~ag,"de  ,S*  if  p-1.2, 

f Ags  d 0 W-  + -±g  S£  if  £ = “ p 

J 2(  1+/0)  2(1 +/9)  1 2 


w 
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